{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nTkZ0i6Yzhym",
    "papermill": {
     "duration": 0.007255,
     "end_time": "2021-01-04T05:10:58.373630",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.366375",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Vaccine RCT Examples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wGGzchhjzhyy",
    "papermill": {
     "duration": 0.006053,
     "end_time": "2021-01-04T05:10:58.386521",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.380468",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "  \n",
    "# Polio RCT\n",
    "\n",
    "One of the earliest randomized experiments were the Polio vaccination trias conducted by the Public Health Service in 1954.  The question was whether Salk vaccine prevented polio.  Children in the study were randomly assigned either a treatment (polio vaccine shot) or a placebo (saline solution shot), without knowing which one they received. The doctors in the study, making the diagnosis, did not know whether a child received a vaccine or not. In other words, the trial was a double-blind, randomized control trial.  The trial had to be large, because the rate at which Polio occured in the population was 50 per 100,000.  The treatment group saw 33 polio cases per 200,745; the control group saw 115 cases per 201,229. The estimated avearage treatment effect is about\n",
    "$$\n",
    "-40\n",
    "$$\n",
    "with the 95% confidence band (based on approximate normality of the two sample means and their differences) is:\n",
    "$$[-52, -28].$$\n",
    "The confidence suggests that the Polio vaccine **caused** the reduction in the risk of polio.\n",
    "\n",
    "The interesting thing here is that we don't need the underlying individual data to evaluate the effectivess of the vaccine. This is because the outcomes are Bernoulli random variabales, and we have enough information to compute the estimate of ATE as well as the confidence intervals.\n",
    "\n",
    "\n",
    "We also compute Vaccine Efficacy metric, which (I googled  [CDC](https://www.cdc.gov/csels/dsepd/ss1978/lesson3/section6.html)) refers to the following measure:\n",
    "$$\n",
    "\\operatorname{VE} = \\frac{\\text{Risk for Unvaccianted - Risk for Vaccinated}}{\\text{Risk for Unvaccianted}}.\n",
    "$$\n",
    "It describes the relative reduction in risk caused by vaccination.\n",
    "\n",
    "\n",
    "It is staighborward to get VE estimate by just plugging-in the numbers, but how do we get the approximate variance estimate. I am too lazy to do calculatios for the delta methods, so I will just use a simulation (a form of approximate bootstrap) to obtain the confidence intervals.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "id": "n0AUYEIczhy2",
    "outputId": "651403dd-4376-4013-dc2b-06ae2656eb14",
    "papermill": {
     "duration": 1.94521,
     "end_time": "2021-01-04T05:11:00.339761",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.394551",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall vaccine efficacy is 0.7124\n",
      "Incidence per 100000 among treated: 16.4388\n",
      "Incidence per 100000 among controlled: 57.1488\n",
      "Estimate of effect on occurances per 100,000: -40.7101\n",
      "Standard error for ATE is: 6.047412320702958\n",
      "95% confidence interval is [-52.5630, -28.8571]\n"
     ]
    }
   ],
   "source": [
    "# Here we calculate the overall effectiveness of the vaccine and construct confidence intervals for it\n",
    "NV =  200745 # number of vaccinated (treated)\n",
    "NU =  201229 # number of unvaccinated (control)\n",
    "RV = 33 / NV # average outcome for vaccinated \n",
    "RU = 115 / NU # average outcome for unvaccinated\n",
    "VE = (RU - RV) / RU # vaccine efficacy\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "\n",
    "print(f\"Incidence per 100000 among treated: {RV * 100000:.4f}\")\n",
    "print(f\"Incidence per 100000 among controlled: {RU * 100000:.4f}\")\n",
    "\n",
    "# treatment effect: estimated reduction in incidence per 100k people\n",
    "effect = 100000 * (RV - RU)\n",
    "\n",
    "print(f\"Estimate of effect on occurances per 100,000: {effect:.4f}\")\n",
    "\n",
    "# variance, standard deviation and confidence interval of ATE using\n",
    "# the fact that outcomes are Bernoulli \n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "var_effect = (100000**2) * (var_rv + var_ru)\n",
    "std_effect = np.sqrt(var_effect)\n",
    "\n",
    "print(f\"Standard error for ATE is: {std_effect}\")\n",
    "\n",
    "ci_effect = [effect - 1.96*std_effect, effect + 1.96*std_effect]\n",
    "\n",
    "print(f\"95% confidence interval is [{ci_effect[0]:.4f}, {ci_effect[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval of VE is [0.5892, 0.8151]\n"
     ]
    }
   ],
   "source": [
    "# we use approximate bootstrap to calculate confidence\n",
    "# intervals for vaccine efficacy via Monte Carlo draws\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV  + np.random.normal(0, 1, B) * np.sqrt(var_rv)\n",
    "RUs = RU  + np.random.normal(0, 1, B) * np.sqrt(var_ru)\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hc5Zn38e8tyeq92Fa1JPfeBDahQyghBUgghIQEEjakZ5Pdlxc24d1NsiRhk03bzW4CCWA6BDBgQosLtnHDyL3JTW6yum3JsmRZI+l5/zhHyjAeSSNpZs6M5v5c11yaOWfmOT+dOTP3nOc0McaglFJKAUQ5HUAppVTo0KKglFKqlxYFpZRSvbQoKKWU6qVFQSmlVC8tCkoppXppUYhwInJGREqdzjEQESkWESMiMfbjt0TkTj+1famI7HV7fFhEPuqPtu32donIFf5qz8dpiog8LiKnRGRjMKc9HJ7vhQo+LQoBJCIr7Q9lnNNZ+mKMSTbGVPq7XRG5S0TW+LvdHsaYjxljnvAhhxGRCQO09Z4xZrI/conIIhF50KP96caYlf5ofxAuAa4BCowxF7qPEJGLRKRVRFI8XyQiW0Tk225F+IzH7TZ/hvR8f/z5XniZ1nnvjTqfFoUAEZFi4FLAAJ8K4HRiAtV2JBjB828ccNgY0+o5whizHqgCPuM+XERmANOA59wGp9s/HHpuLwQytAoBxhi9BeAG/CuwFvg18FePcYuAPwJLgRZgFTDObbwBvgtUAo3AL4Eoe9xddru/AU4CD2IV9weAI0A98CSQZj//NrudVPvxx4BaIMdtWhPccv0v8BZwxp7OWOC3wCmgApjrlvN+4KD9P+wGbraHTwXagS67nSZ7eBzwn8BRoM6eBwl9zL9o+7mNdv5v2Vlj7PErgX+w70+w52Gz/fwX7OGr7de02jluA67A+kK8z54PT/UMc5v2YeBf7P/pFPA4EO82/9d4ZDV2hnsAF9BhT+91t/Y+6jYPfgtU27ffAnH2uJ5s/2y/jzXAl/tZxvKAJfZycAD4qj38bo/5/2Mvr/0BsMJj2C+Axfb9Yvf57cPyngY8amc+jrVcRg/1/fF4L+4FttvPexQYg7WMtgDLgAy3579ov6/NdvvT7eF9vTd5wMtAA3AI+K5bWxcC5cBprOX1105/rwTlu8vpACP1Zn9IvwnMtxfGMW7jFtkL9GX2l8Tv3L9o7A/Ku0AmUATs4+9fgHcBncB3gBggAfiKPb1SIBlYDDzl1t4z9jSzsL6IPuExLfei0GhnjgdW2B+UL2F9ST8IvOv22lvtD1WU/YFuBXLdcnp+ef4W60ssE0gBXgd+3sf8+zpWESq0n/8ufReF54Af2jnigUu8/X/24yvs+fcf9rxPwPsX0U63aa8FHuzn//Kchw96jD/M34vCT4ANwGggB1gH/LtHtp8Ao4AbgDbcvvQ82l2FVcTjgTlYX2xX95XT47WFWMtlkf04Cqsg3WQ/LmZwReFV4GEgyf7fNgJfG+L74/lebMAqBPlYxXIzMNd+/1YA/+b2/K9gLVs9xXerx+fuQbfHUcAmrB9wsVifn0rgOnv8euCL9v1kYKHT3yvBuDkeYCTesPpzXUC2/bgC+L7b+EXA826Pk7F+1RXajw1wvdv4bwLL7ft3AUc9prcc+Kbb48n29Hu+QNOxfp3vAB72eK3nF9qf3MZ9B9jj9ngm9q/+Pv7vrcCNbjndC51gFY3xbsMuAg710dYK4Otuj6+l76LwJPAIVv+5ZzvevnQ6sH/5uw3z/CJyn/YNwEFv/1cf87C/onAQuMFt3HVY3Tw9Oc7i9kWM9SV43pcR1pd6F5DiNuznwKK+cnppYxnwA/v+NVg/CEbZj4vt/6vJ4zbVSztjgHO4rfUBt2P/gBjC++P5XnzB7fHLwB88ltFX+/j/0u32e9aaP/TeAAs4/7P0L8Dj9v3VwI+xP8eRctNtCoFxJ/A3Y0yj/fhZe5i7Yz13jDFnsLoA8ryNx+oW6msc9rgjHs+PwfqwYoxpwlqtngH8aoDsdW73z3p5nNzzQES+JCJbRaRJRJrs9rP7aDcHSAQ2uT3/bXu4N3mcPw/68n+xis5Ge0+fr/TzXIAGY0z7AM/pb/4Ph7f3yr3tE8aYTrfHbbjNc492ThpjWjzayh9Eliew1gIBvgg8a4xxeTwn2xiT7nbb46WdcVhrNjVu7+3DWGsMMPj3x5NPy6SIRIvIQyJyUEROYxUU6HuZHAfk9WS2c/8A+3OD1Q03CagQkQ9E5BODzB2WRupGNseISALwWSBaRGrtwXFAuojMNsZss4cVur0mGaubotqtqUJgl32/yGOc8ZhsNdYC3qMIqxuizm5/DtZq9XPAfwHXD+mfcyMi44A/AVcD640xXSKyFevD7y1jI9YHeLox5rgPk6jBbR5h/U9eGWNqga/auS4BlonIamPMgb5e4sP0PafdM/9bsYob9vTGDrLtnveqr/fWV9VApoikuBWGIqz+fF8tBv5XRK4EPo31K30ojmGtKWR7FDRgSO/PUH0euBH4KFZBSMPaJtTXMnkMa011orfGjDH7gdtFJApr/rwkIlnGy8b7kUTXFPzvJqzV+mlY/bxzsDa8vsfff5UB3CAil4hILPDvwPvGGPdfp/eKSIaIFAL/CPS318dzwPdFpMQuMD/D2pjXKSLxwNNYv4C+DOSLyDf98H8mYX3IGgBE5MtYawo96oAC+//DGNONVUR+IyKj7dfki8h1fbT/F+C7IlIgIhlYG7W9EpFbRaTAfnjKztXllmMox2F8y552Jta865n/24DpIjLHnrc/8njdQNN7DnhARHJEJBurP/vpwYazl5V1wM9FJF5EZmH9sn1mEG20Ai9hbUg/YowpH2wOu50a4G/Ar0QkVUSiRGS8iFwOAXt/vEnBKk4nsAr3zzzGe05rI3BaRO4TkQR7TWOGiFxg575DRHLsZbfJfk0XI5wWBf+7E6tP8qgxprbnBvwe+ILbLpDPAv+G1W00H/iCRzuvYW0E2wq8gbXXRV8ew9qLZjXWhuF2rL5WsPqZq4wxfzDGnAPuAB4UEa+/jnxljNmN1RW1HuvDNhNrg2yPFVi/hmtFpKcb7T6sDeIb7NX7ZVjbP7z5E/AO1pfwZqxftX25AHhfRM5gbcj+R2PMIXvcj4An7O6Bzw7iX3wW64uu0r49CGCM2Ye1IXgZsB/wPBbjUWCaPb1XvbT7INYeLduxtvFs7ml7CG7H6vuvBl7B2uC6dJBtPIG15vJkH+ObPI5T+Kc+nvclrI21PXtsvQTk2uMC8f548yRWF9pxO8cGj/Efem+MMV3AJ7F+uB3CWpv9M9YaBlhr1Lvs3L8DPudDt2PYE3uDigoiEVmE9UX9QB/jDTAxAKvXSinVL11TUEop1UuLglJKqV7afaSUUqqXrikopZTqFRbHKWRnZ5vi4mKnYyilVFjZtGlTozGmrwNEvQqLolBcXEx5+ZB2oVZKqYglIv2dCcAr7T5SSinVS4uCUkqpXloUlFJK9dKioJRSqpcWBaWUUr20KCillOqlRUEppVQvLQpKKaV6hcXBa0qNFO2uLp7ZcIQPjpzibEcXl03K4dNz88lIinU6mlKArikoFTRbjzVxxS9X8rfddYzPSWZOYTqr99VzzW9W8W5FvdPxlAJ0TUGpoFh3sJFvPrOZuy8poWxcZu/whaVZ7K45zb0vbeP+66dwS1lhP60oFXhaFJQKsIMNZ/jmM5v5zpUTmJaXdt74abmp3P+xqfz8zT3EjYrmk7PzHEiplEW7j5QKoNZznXz1yXJunV/gtSD0yE9P4N7rJvPAqzvZebw5iAmV+jAtCkoF0H+8XUFhRgJXTRkz4HPHZSVx50XFfPXJcprbXEFIp9T5tCgoFSDlh0/yxvYa7lhY7PNrLhqfxZzCdO5bvB29KqJyghYFpQKgs6ub+17ezh0Lx5EcN7hNd5+7oIg9Nad5eVNVgNIp1TctCkoFwHMbj5IUG8OCksyBn+whNiaKb1w+np++uYfa5vYApFOqb1oUlPKz0+0ufrNsP7cvKEJEhtTGuKwkrp46hvte1m4kFVxaFJTys4dXVTK7II3irKRhtXPj7DyOnmxlybZqPyVTamBaFJTyo1OtHTy1/jA3zckfdlsx0VF85eISfvL6bpraOoYfTikfaFFQyo8eXn2QBSWZjE6N90t7E0anUFacwc/frPBLe0oNRIuCUn7SfNbFM+8f5ZOzh7+W4O7W+YUs3VPH9qomv7arlDdaFJTyk6c3HGFeUQY5KXF+bTcpLoZb5xfwwCs76e7Wjc4qsAJWFESkUETeFZE9IrJLRP7RHp4pIktFZL/9NyNQGZQKlnZXF4+vPcQNM3MD0v5lk3I46+ri9e260VkFViDXFDqBfzbGTAUWAt8SkWnA/cByY8xEYLn9WKmwtmRbNUWZiRRlJgak/SgRPndhEQ+9VUG7qysg01AKAlgUjDE1xpjN9v0WYA+QD9wIPGE/7QngpkBlUCoYjDE8vvYQ10wb+PxGwzEtN5WCjASeXH84oNNRkS0o2xREpBiYC7wPjDHG1IBVOIDRfbzmHhEpF5HyhoaGYMRUaki2HGuiuc3FrIL0gE/r1vmF/GHlQVrPdQZ8WioyBbwoiEgy8DLwPWPMaV9fZ4x5xBhTZowpy8nJCVxApYZp0drDXDV1NFFDPHp5MAozE5mam8rTG44EfFoqMgW0KIjIKKyC8IwxZrE9uE5Ecu3xuYBeh1CFreY2Fysq6rlsYvB+uNw4J5+HV1dytkO3LSj/C+TeRwI8CuwxxvzabdQS4E77/p3Aa4HKoFSgvb69mlkFaaTEjwraNIsyE5kwOpmXNh0L2jRV5AjkmsLFwBeBq0Rkq327AXgIuEZE9gPX2I+VCkvPf3CUSydmB326108fy5/fO6THLSi/C9g1mo0xa4C+OlmvDtR0lQqW/XUt1DS1MzM/8BuYPU0Zm0JMtLByX71PV3VTyld6RLNSQ7R483EunpBNdFTgNzB7EhGus9cWlPInLQpKDYExhiXbqllYmuVYhgUlWeyqPs2xk22OZVAjjxYFpYZgx/FmDIbirMAcweyL2JgoLp6QxXMbjzqWQY08WhSUGoIl26pZUJI15Cur+csVk0bzYnkVnV3djuZQI4cWBaUGyRjDG9trHO066lGYmUhmciyr9+tR/8o/tCgoNUi7qk8jAoUZCU5HAeAj47N4edNxp2OoEUKLglKD9M6uWuYVZTjeddRjYUkWK/fV6/mQlF9oUVBqkN7ZVcv8otC5DEhqwiimjEll6e46p6OoEUCLglKDcOxkG/WnzzFpTIrTUT5k4fgsFm+ucjqGGgG0KCg1CMv31DFvXAZRDhyw1p/5RRmUHzlF81mX01FUmNOioNQgLNtTz6yCNKdjnCchNprpeamsqNAuJDU8WhSU8lG7q4tNR08xIy/0igLA/HGZ/HVbjdMxVJjToqCUjzZUnqA0O4mkuICdR3JY5o/LYH3lCd0LSQ2LFgWlfPRuRT0z8kNzLQEgOS6GyWNSWLlXD2RTQ6dFQSkfrdzbwOwgXId5OOYWpfPOrlqnY6gwpkVBKR8cO9lGc7uLcQ6eAM8X84oyWLm3HpeeC0kNkRYFpXyw7mAjM/PTiAqRo5j7kpUcx+jUeMoPn3I6igpTWhSU8sGqfQ1MzU11OoZP5mkXkhoGLQpKDaC727D+4AlmhvBGZnfzijJYursOY/T6zWrwtCgoNYCK2hYSYqPJTo5zOopPijIT6ejq5mDDGaejqDCkRUGpAaw90BiyB6x5IyLMK0xn2Z56p6OoMKRFQakBvLc/fLYn9JhdmK5nTVVDokVBqX50dnWz6egppoVZUZiel8aemtM0tXU4HUWFGS0KSvVjx/FmclLiSE0Y5XSUQYmNiWJ6Xiqr9unRzWpwtCgo1Y/1B08wdWx4rSX0mF2gXUhq8LQoKNWPNQcaw257Qo85hems3tdApx7drAZBi4JSfejo7GbrsaawLQpZyXFkp8Sx5ViT01FUGNGioFQfdhxvYmxaPMkheqpsX8wuSGf5Hu1CUr7ToqBUH8J5e0KPuYXpLNutxyso32lRUKoPaw+cYEpuitMxhmV8TjKNZ85xvOms01FUmNCioJQXHZ3dbK1qCvs1hagoYXZhOisqdG1B+UaLglJebK9qIj89IWQvvTkYswvSWKpnTVU+0qKglBfrK08wZWx4dx31mFWQTvmRU5zt6HI6igoDWhSU8mLN/kamhHnXUY+kuBhKc5JYd7DR6SgqDGhRUMpDu6uL7cebmRrmG5ndWbum6nYFNTAtCkp52Hz0FEUZiSTGhv/2hB5z7I3NeuEdNZCAFQUReUxE6kVkp9uwH4nIcRHZat9uCNT0lRqqdQdOjKi1BID89ATAsK9OL7yj+hfINYVFwPVehv/GGDPHvr0ZwOkrNSRrDjQyLYwuquMLkZ5dU/XoZtW/gBUFY8xq4GSg2lcqEM6c62RvbQuTxiQ7HcXv5uiFd5QPnNim8G0R2W53L2X09SQRuUdEykWkvKFBzwmvguP9yhNMGJ1MXEy001H8blpuGhW1LTS3uZyOokJYsIvCH4DxwBygBvhVX080xjxijCkzxpTl5OQEK5+KcKv2NTAjf2TsiuopNiaKabmpvHdAf2SpvgW1KBhj6owxXcaYbuBPwIXBnL5SA1m1t4GZ+elOxwiYmQVpetZU1a+gFgURyXV7eDOws6/nKhVsx0620dzuYlxWotNRAmZOQTqr9jbQ3a27pirvArYjtog8B1wBZItIFfBvwBUiMgcwwGHga4GavlKDteZAI7Py04gScTpKwIxOjScpLoad1c3MKhi5a0Rq6AJWFIwxt3sZ/GigpqfUcK2oqGdG/sjaFdWb2YXpvFtRr0VBeaVHNCuFdWqL9QdPMLtw5H9RzirQU2mrvmlRUArrrKhFmYmkxo9yOkrATR6Twv76MzS1dTgdRYUgLQpKAX/bVcvcopG/lgBuu6bu17OmqvNpUVARzxjDst31zCvq81jKEWdmfhrvaheS8kKLgop426uaiY2JIi89wekoQTO7MJ1V+xr0rKnqPFoUVMRbsq2aBSWZTscIqjGp8cTFRLGnpsXpKCrEaFFQEa272/D6tmoWlmY5HSXoZhaksWqfdiGpD9OioCJa+ZFTJMZGU5g5co9i7suM/DRW7tXzIKkP06KgItqrW46zIALXEgCm56axvaqZto5Op6OoEKJFQUWsdlcXb+yo4eLxkVkUEmKjGT86iQ2VJ5yOokKIFgUVsd7ZVUtpdhI5KfFOR3HMjDztQlIfpkVBRaxn3j/KZZMi+1odswqsXVOV6qFFQUWkw42t7K9roWxc5Byw5s24rESa2lxUnWpzOooKEVoUVERatO4Ql0/KISY6sj8CUSLMKkhjjZ7yQtl8+kSIyMsi8nERiexPkBoRWtpdvLz5OB+dOsbpKCFhel6qbldQvXz9kv8D8Hlgv4g8JCJTAphJqYB6sfwYM/PTyEqOczpKSJiZn866g4106dXYFD4WBWPMMmPMF4B5WFdMWyoi60TkyyIy8s81rEaMrm7Do2sOc930sU5HCRmZSbFkJMWy43iz01FUCPC5O0hEsoC7gH8AtgC/wyoSSwOSTKkAeHtnLakJMUwak+J0lJAyIy+NVXv1lBfK920Ki4H3gETgk8aYTxljXjDGfAdIDmRApfzFGMMfVh7gYzNynY4Scmbmp7FSd01V+L6m8GdjzDRjzM+NMTUAIhIHYIwpC1g6pfxo46GTnGpzMT/Cd0P1ZmpuKhU1LZxudzkdRTnM16LwoJdh6/0ZRKlA++Oqg1w3fSxRIk5HCTmxMVFMHpvC+oN6yotI129REJGxIjIfSBCRuSIyz75dgdWVpFRYqGw4w5ajTVw2KdvpKCFrRn6qXo1NETPA+OuwNi4XAL92G94C/CBAmZTyuz+9V8nVU0cTFxPtdJSQNSs/nd8u24cxBtG1qYjVb1EwxjwBPCEinzHGvBykTEr5VXObi9e31fDLW2Y5HSWkFWQk4OrqprKxlfE5uv9IpOq3KIjIHcaYp4FiEfknz/HGmF97eZlSIeXFTceYW5ROemKs01FCmohY127e26BFIYINtKE5yf6bDKR4uSkV0rq7DYvWHdZTWvhoRn4aK3S7QkQbqPvoYfvvj4MTRyn/eu9AI7HRUUwcrb98fTEjL42HV1fS7uoifpRuf4lEvh689gsRSRWRUSKyXEQaReSOQIdTarieff8IV0werRtOfZQUF0NpdpLumhrBfD1O4VpjzGngE0AVMAm4N2CplPKDk60drD1wgo9E6OU2h2pWQRrL9tQ5HUM5xNei0HPSuxuA54wxJwOURym/eXXLceYVpZMUN9Ce18rd3MIMVlTUY4yeNTUS+VoUXheRCqAMWC4iOUB74GIpNXzPf3CUSydG9uU2h6IgI4GubsP++jNOR1EO8PXU2fcDFwFlxhgX0ArcGMhgSg3H3toWTrV2MC0v1ekoYUdEmFuYznLtQopIg7mS2lTgNhH5EnALcG1gIik1fK9tPc7C0iw9z9EQzS5MZ+luLQqRyNe9j54C/hO4BLjAvunZUVVIMsbw2tZqLhqv5zkaqul5aeytbeFka4fTUVSQ+boFrgyYZnTLkwoDW441IQLFWXrOxqGKjYliZoHVhXRrWaHTcVQQ+dp9tBPQ6xeqsPD6tmoWlGTqsQnDNLcwnbd31TodQwWZr0UhG9gtIu+IyJKeW38vEJHHRKReRHa6DcsUkaUist/+q1c7UX5ljOGtHbUsKNFjE4ZrblE6GypP0O7qcjqKCiJfu49+NIS2FwG/B550G3Y/sNwY85CI3G8/vm8IbSvl1faqZmKihYKMBKejhL2U+FGUZifz3v5Grpmm546KFL7ukroKOAyMsu9/AGwe4DWrAc+D3G4EnrDvPwHcNJiwSg3kjR01XDAuQ7uO/GT+uAze2F7tdAwVRL7uffRV4CXgYXtQPvDqEKY3pucaz/bf0f1M8x4RKReR8oYGvaC4Gpgxhjd31HCBdh35zQXFmayoqKejs9vpKCpIfN2m8C3gYuA0gDFmP/18ofuDMeYRY0yZMaYsJ0ePSlUD21d3BldXt+515EeZSbHkpSew7mCj01FUkPhaFM4ZY3p3WBaRGGAou6fWiUiu3UYuoCduV37zzq5a5mvXkd+VFWfwxvYap2OoIPG1KKwSkR8ACSJyDfAi8PoQprcEuNO+fyfw2hDaUMqrt3fWML9Id2jztwuLs/jb7jpcXdqFFAl8LQr3Aw3ADuBrwJvAA/29QESeA9YDk0WkSkTuBh4CrhGR/cA19mOlhq266SxVTWeZPFbPdeRvOSlx5KbFs+aAdiFFAp92STXGdIvIq8CrxhiftvoaY27vY9TVvoZTyldLd9cxtzCD6CjtOgqEC0syeW3Lca6cHNBNiSoE9LumIJYfiUgjUAHsFZEGEfnX4MRTyjdv7axlnnYdBcyCkiyWV9TrgWwRYKDuo+9h7XV0gTEmyxiTCSwALhaR7wc8nVI+aD7rYntVE7MK0pyOMmJlJsVSnJXEyr26e/hIN1BR+BJwuzHmUM8AY0wlcIc9TinHrdxbz7TcVL3QfIAtLM1i8eYqp2OoABuoKIwyxpy3dcnerjDKy/OVCrp3dtYypyjd6Rgj3oKSTNYeaKT5rMvpKCqABioK/Z1MXU+0rhzX0dnNewcadVfUIEiKi2FGfhpv79RjFkaygYrCbBE57eXWAswMRkCl+rOh8gT56QmkJ8Y6HSUiXDQ+i5c2aRfSSNZvUTDGRBtjUr3cUowx2n2kHPfWzhrmj9O1hGCZV5TB3toWjjeddTqKCpDBXKNZqZDS3W1YurtOi0IQjYqOYkFpFq9u0bWFkUqLggpb26qaSIyNITdNr50QTBePz+alTcfRq/OOTFoUVNh6Z1ct83Svo6CbNCaZdlcXO443Ox1FBYAWBRWWrGsn1FJWnOl0lIgjIlw8PouXdYPziKRFQYWlfXVnaHd1UZqd5HSUiHTxhByWbKvWM6eOQFoUVFh6c0cNFxRn6rUTHDI2LZ6xafGs3qenvRhptCiosPTmjhrKinWvIyddVJrNi+XahTTSaFFQYedQYyuNZ84xaXSK01Ei2kWlWby3v4HT7Xrai5FEi4IKO29sr2ZBSSZReu0ERyXH26e92FHrdBTlR1oUVNhZsq2aC0qynI6hsE97oWdOHVG0KKiwUtlwhsYzHUwZo11HoWBuYQYVNaep1tNejBhaFFRYeWNHDRcUZ2jXUYiIjYmyLtW59bjTUZSfaFFQYWXJ1moWaNdRSLmoNItXtmhRGCm0KKiwsa+uhaa2DiaP1a6jUDJlbConWzvYV9fidBTlB1oUVNhYsrWaBaVZROkBayElKkpYWJrFa7q2MCJoUVBhwRjDq1uPs7BUu45C0UfGZ/Pq1mo9c+oIoEVBhYXtVc10G6PnOgpRxVmJiMCWY01OR1HDpEVBhYXFm6u4qDRLz3UUokS0C2mk0KKgQl5nVzevb6/h4vHZTkdR/fhIaRZ/3V5DV7d2IYUzLQoq5K09eILs5Fhy0/UKa6EsNz2BjKRY1h884XQUNQxaFFTIW7y5Sjcwh4mFpZm8ql1IYU2Lggpprec6Wb6nnou0KISFi0qzeWd3Le2uLqejqCHSoqBC2ju7apk8NoX0xFinoygfZCbFUpKVxMq9evGdcKVFQYW0lzZVcfF4XUsIJwtKs3hli545NVxpUVAhq7a5nR3Hm5k/LtPpKGoQLizJZM3+RprP6sV3wpEWBRWyFm+uYkFJJrExupiGk+Q46+I7b+2ocTqKGgL9tKmQZIzhhQ+OcenEHKejqCG4eHw2L27SLqRwpEVBhaTNR0/RZQwTRyc7HUUNwZyidPbXt3DsZJvTUdQgaVFQIen5jce4dGK2ntYiTI2KjuKikixe2azHLIQbR4qCiBwWkR0islVEyp3IoELXmXOdvL2zVruOwtylk3J4ofwY3Xrai7Di5JrClcaYOcaYMgczqBD02tbjTM9PJUOPTQhrpdlJjIoWNh4+6XQUNQjafaRCijGGp9Yf4YpJo52OooZJRLh0Yg7PbzzqdBQ1CE4VBQP8TUQ2icg93p4gIveISLmIlDc06F+IgzAAAA+TSURBVNGRkWJbVTNNbR3MLEhzOoryg0smZLNsT70esxBGnCoKFxtj5gEfA74lIpd5PsEY84gxpswYU5aTo33LkeLRNZVcPXWMXnJzhEhNGMWsgjQWb9bdU8OFI0XBGFNt/60HXgEudCKHCi31p9t5t6KBKyZr19FIcvWU0Ty1/oheqjNMBL0oiEiSiKT03AeuBXYGO4cKPU+tP8LF47NIjotxOoryo6m5qXR2d/PB4VNOR1E+cGJNYQywRkS2ARuBN4wxbzuQQ4WQ1nOdPP3+Ea6dPtbpKMrPRISrpozhsTWHnI6ifBD0n2TGmEpgdrCnq0LbM+8fYWpuKnl6dbUR6bKJOXzvhS0cbzpLvr7HIU13SVWOa3d18cjqSj41O8/pKCpAEmKjuWxSDovWHnY6ihqAFgXluKc3HKEkO4lxWUlOR1EBdO20Mfyl/BhnznU6HUX1Q4uCclRzm4v/efcAt84vdDqKCrCclHhm5qfxzIYjTkdR/dCioBz1+3cPMH9cBoWZiU5HUUHwiVm5/Pm9Q3oN5xCmRUE55kD9Gf5Sfoyb5xY4HUUFybisJIqyEnlZr7UQsrQoKEcYY/jB4u3cOCePzCQ98V0k+dTsPP57xQHOderaQijSoqAc8WL5MU60dnDtND0uIdJMGpNCXnoCL2w85nQU5YUWBRV0Vafa+NmbFXz10lKio/QcR5Ho0/Py+e8VB3TbQgjSoqCCqqvb8P0XtnLDzLG6C2oEG5+TzPjRSTyqRzmHHC0KKqj+a/l+zrq6+PhMPVAt0n22rJBHVldy4sw5p6MoN1oUVNCsO9DIUxuO8I3LJxCl3UYRLzctgY+Mz+JXf9vndBTlRouCCoqa5rN857ktfP3y8bq3ker16bkFvLWzhl3VzU5HUTYtCirgznV28bWnNnHt9DHMzNcrqqm/S46P4Zb5BfzwlZ16vYUQoUVBBZQxhh8u3klibDSfnKXbEdT5rpg0mtZznTz/ge6iGgq0KKiAemLdYcqPnORrl41H9BKbyouoKOHuS0r4xdsV1J9udzpOxNOioAJm7YFGfrd8P9//6CTiR0U7HUeFsHFZSVwxebR2I4UALQoqIA41tvLtZzfz7SsnMDo13uk4KgzcNCefvXUtLNlW7XSUiKZFQfldc5uLux7fyC3zC5iWpxuWlW9iY6L42mWl/GjJLmqbtRvJKVoUlF91dHZzz1PlzMhP46opY5yOo8JMaU4yH506hu+9sJXubu1GcoIWBeU3xhjufWkbAJ+/oMjhNCpc3Tgnn9NnXTy8+qDTUSKSFgXlNz9/q4KKmha+ccV4PWJZDVl0lPD1y8fzyOpKNh055XSciKNFQfnFH1ce5K0dNfzztZOIi9E9jdTw5KTEcfclpXzrmc2cbO1wOk5E0aKghu1Pqyt5Yv1h7rt+Cinxo5yOo0aI+eMyuLAkk28+s5nOrm6n40QMLQpqyIwx/H7Ffh5be4gf3DCVrOQ4pyOpEea2skLaXV08+MYep6NEDC0Kakg6u7r519d28eKmKh74+DSytSCoAIiKEr515QSW7alj0Vq99kIwxDgdQIWfhpZzfPvZzXR0dvPAx6eRHKeLkQqc5LgY7r12Mj/5625yUuL5+KxcpyONaPppVoOydHcd97+8ncsn5fCZeQW6l5EKitGp8dx73WQeeHUH0VFw/QwtDIGiRUH55NjJNn7y+m5215zm21dOYEpuqtORVIQZl5XEvddN4V8W7+BUm4vbL9RjYQJBi4LqV+OZc/zPuwd4eVMV188Yy89unklsjG6KUs4oyU7i/31iGr98Zy/761q4/2NTdXn0M52byqtTrR089NYervrPldQ0neU/PjOLm+cW6AdQOS43LYEff2o6O4438+n/Xcve2hanI40ouqagPqS5zcUj71Xy1PrDLCjJ5Kc3z9Q9i1TISYkfxfc/OokVFfV89uH1fGZePt+9eiLpiXqp1+GScDh3eVlZmSkvL3c6xojW0u7i8TWHeXTtIeaPS+emOfnkpOgpr1Xoa2rrYPGW42w8dJIvLRzHly8p0euA20RkkzGmbFCv0aIQ2U63u3hi3WEeW3OIGflp3Dw3n9y0BKdjKTVodafb+ev2at4/dJJPzsrjK5eUMGF0stOxHKVFQfns6Ik2nlh/mJc2VTG7MJ1PzcojP0OLgQp/TW0dLNtTx/KKeqblpnLXR4q5aspoYqIjb3uYFgXVr/qWdpburuO1rdXsrW3h0onZXDttjHYTqRGpo7Ob9w+dYEVFPSdbO7jtgkI+W1ZIYWai09GCRouC6mWM4djJs2yrauKDwydZf/AENc3tzClMp6w4g3lFGYyKwF9OKjIdOdHKqn0NrDt4gpLsJG6ak8e108eSlz6y147DpiiIyPXA74Bo4M/GmIf6e74Whf6dau1gX10L++rPUFFzmj01p9lXd4b4UVGMz0mmNCeJKWNTGZ+TTLQegawiWGdXN9urmvng8Ek2Hz1FTkocF0/Ipqw4kxl5qYzLShpRn5GwKAoiEg3sA64BqoAPgNuNMbv7eo0WBXB1dVPddJYjJ9qobDjDvroz7K9v4UD9GTo6uynMTKQgI4HctASKMhMpykwkNUFPY61UX7q7DYdOtLK7+jQHG85w5GQbp1o7KMhIoCAjkfz0BMamxZOTEkdOchzZKXHkpMSRnRwbNtcMGUpRcOI4hQuBA8aYSgAReR64EeizKAxH81kXdafbMQYMxvrrdr+r29BlDN3dxrpvP+7qtsZ/iIBYmen5LdFlDF1dBldXN+2dXZzt6Kato5Mz5zppae+kpd1FS3snrec6Oevqot3Vjaurm65ugwhEiRAlQky09VfEyufq6qato4uWdhenWl10GUNcTBRj0+LJTYtnbGoCC0qy+PS8AjISY/H22+b0WVcgZqlSI0ZOchyXT8rh8kk5AJzr7KL29DkaWs5xovUcW481cfqsi2b7dqqtA1eX9cWQljCKrKRYspJjyUyKJSMxlrSEUSTFxZAYG038qGjiYqKIjYkiNjqKmOgoYqKFUVHW35goISpKiBYhOsr67AtCVJT1F0AESrOTgrqR3ImikA8cc3tcBSzwfJKI3APcA1BUNPRznHzj6U2sO3hiyK8PFSnxMUSL0NzmornNRUWNHsWpVLAlxcZg4EOForKxNaDT/MUts/hsWWFAp+HOiaLg7UfteX1YxphHgEfA6j4a6sSe/erCob5UKaUijhO7n1QB7mWvAKh2IIdSSikPThSFD4CJIlIiIrHA54AlDuRQSinlIejdR8aYThH5NvAO1i6pjxljdgU7h1JKqfM5cpZUY8ybwJtOTFsppVTf9JBWpZRSvbQoKKWU6qVFQSmlVC8tCkoppXqFxVlSRaQBOOIxOBtodCDOcGnu4NLcwaW5g2ug3OOMMTmDaTAsioI3IlI+2BM9hQLNHVyaO7g0d3AFIrd2HymllOqlRUEppVSvcC4KjzgdYIg0d3Bp7uDS3MHl99xhu01BKaWU/4XzmoJSSik/06KglFKqV0gUBRG5XkT2isgBEbnfy/giEXlXRLaIyHYRucEeXiwiZ0Vkq337o9tr5ovIDrvN/xIRv1+Nexi5v+CWeauIdIvIHHvcSrvNnnGjHcg9TkSW25lXikiB27g7RWS/fbvTbXgozG+vuUVkjoisF5Fd9rjb3F6zSEQOuc3vOaGS2x7X5ZZtidvwEhF5334fXrBPQx8SuUXkSo/lu11EbrLHBWN+PyYi9SKys4/xYi+jB+zs89zGObl8Dym335dvY4yjN6zTZx8ESoFYYBswzeM5jwDfsO9PAw7b94uBnX20uxG4COtKb28BHwuV3B7PmQlUuj1eCZQ5PL9fBO60718FPGXfzwQq7b8Z9v2MEJrffeWeBEy07+cBNUC6/XgRcEsozm/78Zk+2v0L8Dn7/h97lrNQye32nEzgJJAYjPltT+MyYF4/3w032MuoAAuB951evoeZ26/LdyisKVwIHDDGVBpjOoDngRs9nmOAVPt+GgNcqU1EcoFUY8x6Y82ZJ4Gb/Bvbb7lvB57zc7b++JJ7GrDcvv+u2/jrgKXGmJPGmFPAUuD6EJrfXnMbY/YZY/bb96uBemBQR3kOw3Dmt1f2r9SrgJfsQU8QQvPbwy3AW8aYNj/n65MxZjVWIerLjcCTxrIBSLeXYSeX7yHn9vfyHQpFIR845va4yh7m7kfAHSJShXUdhu+4jSsRq3tmlYhc6tZm1QBtDtdwc/e4jfOLwuP2qt7/C8Bqqi+5twGfse/fDKSISFY/rw2V+d1X7l4iciHWL9+DboN/aq92/0ZE4vwbe9i540WkXEQ29HTBAFlAkzGms582nc7d43Ocv3wHcn77or/l2Knl2xcDvif+WL5DoSh4+9Lz3E/2dmCRMaYAaxXqKRGJwlpNKjLGzAX+CXhWRFJ9bHO4hpPbakBkAdBmjHHvQ/yCMWYmcKl9+6J/Y/uU+/8Al4vIFuBy4DjQ2c9rQ2V+95XbasD6xfcU8GVjTLc9+F+AKcAFWN0G94VY7iJjncbg88BvRWS8j20Ol7/m90ysqyz2CPT89sVgl+NgzG9f9JvDX8t3KBSFKqDQ7XEB53ez3I3Vh4oxZj0QD2QbY84ZY07YwzdhVcdJdpsFbq/31qZjud3Gn/cryhhz3P7bAjyLtRrvTwPmNsZUG2M+bRfbH9rDmvt5bUjM735yY/9YeAN4wF717nlNjb06fg54nNCa3z3dARhjKrG2N83FOgFauojE9NWm07ltnwVeMca43F4T6Pnti/6WY6eWb1/0+Z74dfkezAaIQNywLglaCZTw9w1a0z2e8xZwl31/qj0jBKvfLNoeXor1SyXTfvwB1saYng1DN4RKbvtxlP0ml3q0mW3fH4XVZ/x1B3JnA1H2/Z8CP7HvZwKHsDbCZdj3Q2l+95U7Fqvv+3te2s21/wrwW+ChEMqdAcS5PWc/9sZerI287huavxkqud3GbwCuDOb8dptOMX1vsP04H95gu9Hp5XuYuf26fPv9nxrijLgB2If1S/+H9rCfAJ+y708D1toL5lbgWnv4Z4Bd9vDNwCfd2iwDdtpt/h77yzgUctvjrgA2eLSXBGwCttv/1++wi16Qc9+C9QW0D/gz9heTPe4rwAH79uUQm99ecwN3AC77Pei5zbHHrQB22NmfBpJDKPdH7Gzb7L93u7VZirVHzAGsAhEXKrntccVYP9KiPNoMxvx+Dqtr2YX1w+tu4OvYP7CwviD/x/6/duC2t5/Dy/eQcvt7+dbTXCillOoVCtsUlFJKhQgtCkoppXppUVBKKdVLi4JSSqleWhSUUkr10qKglFKqlxYFpZRSvf4/y+n09HIxDhcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "79c7e3d0-c299-4dcb-8224-4455121ee9b0",
    "_uuid": "d629ff2d2480ee46fbb7e2d37f6b5fab8052498a",
    "id": "In0Vkgi1zhy4",
    "papermill": {
     "duration": 0.007763,
     "end_time": "2021-01-04T05:11:00.355927",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.348164",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Pfizer/BNTX Covid-19 RCT\n",
    "\n",
    "Here is a link to the FDA [briefing](https://www.fda.gov/media/144245/download) and an interesting [discussion](\n",
    "https://garycornell.com/2020/12/09/statistics-in-the-pfizer-data-how-good-is-the-vaccine/?fbclid=IwAR282lS0Vl3tWmicQDDhIJAQCMO8NIsCXyWbUWwTtPuKcnuJ2v0VWXRDQac), as well as data.\n",
    "\n",
    "Pfizer/BNTX is the first vaccine approved for emergency use to reduce the risk of Covid-19 decease. Volunteers were randomly assigned to receive either a treatment (2-dose vaccination) or a placebo, without knowing which they recieved. The doctors making the diagnoses did not know now whether a given volunteer received a vaccination or not. The results of the study are given in the following table:\n",
    "\n",
    "![pfizer/biontech result table](https://lh6.googleusercontent.com/oiO6gYom1UZyrOhgpFx2iq8ike979u3805JHiVygP-Efh1Yaz2ttyPcgWKlT1AqHDM4v46th3EPIkOvRLyXA0fNUloPL-mL9eOFmSAzfbNOHyCZSQ0DyzMhcFUtQuZ520R5Qd2lj)\n",
    "\n",
    "Here we see both the overall effects and the effects by age group. The confidence intervals for the averal ATE are tight and suggest high effectivness of the vaccine. The confidence intervals for the age group 65-75 are much wider.  We could group 65-75 and >75 groups to evaluate the effectiveness of the vaccine and also narrow down the width of the confidence band. \n",
    "\n",
    "In this case, the reported results are for vaccine effectiveness. We use the same approach as above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wLulUAGdzhy5",
    "papermill": {
     "duration": 0.007593,
     "end_time": "2021-01-04T05:11:00.371394",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.363801",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the code cell below  we calculate the overall effectiveness of the vaccie and construct confidence intervals for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "fmNIIhxuzhy9",
    "outputId": "34cc5e18-729b-4edc-f489-96a392068c1a",
    "papermill": {
     "duration": 0.497816,
     "end_time": "2021-01-04T05:11:00.876971",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.379155",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall vaccine efficacy is 0.9462\n",
      "Incidence per 100000 among treated: 45.0789\n",
      "Incidence per 100000 among controlled: 837.7950\n",
      "Estimate of effect on occurances per 100,000: -792.7161\n",
      "Standard error for ATE is: 65.91015718590171\n",
      "95% confidence interval is [-921.9000, -663.5322]\n"
     ]
    }
   ],
   "source": [
    "NV =  19965\n",
    "NU =  20172\n",
    "RV = 9 / NV\n",
    "RU = 169 / NU\n",
    "VE = (RU - RV) / RU\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "\n",
    "print(f\"Incidence per 100000 among treated: {RV * 100000:.4f}\")\n",
    "print(f\"Incidence per 100000 among controlled: {RU * 100000:.4f}\")\n",
    "\n",
    "# treatment effect: estimated reduction in incidence per 100k people\n",
    "effect = 100000 * (RV - RU)\n",
    "\n",
    "print(f\"Estimate of effect on occurances per 100,000: {effect:.4f}\")\n",
    "\n",
    "# variance, standard deviation and confidence interval of ATE using\n",
    "# the fact that outcomes are Bernoulli \n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "var_effect = (100000**2) * (var_rv + var_ru)\n",
    "std_effect = np.sqrt(var_effect)\n",
    "\n",
    "print(f\"Standard error for ATE is: {std_effect}\")\n",
    "\n",
    "ci_effect = [effect - 1.96*std_effect, effect + 1.96*std_effect]\n",
    "\n",
    "print(f\"95% confidence interval is [{ci_effect[0]:.4f}, {ci_effect[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval of VE is [0.9086, 0.9809]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEICAYAAAC55kg0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gc5bn38e+9WvXerF7ce7epJmB6CRBOSAgJCZBCEkjPm5P6vgkJOScnJ42TnBQIwZTgkFASCMUYG2PccG9ylW3Jli3JkizZlmz15/1jRmFZtPJK2t3Z1d6f69pLuzO7z/w0O7v3zDOzM2KMQSmlVHRyOR1AKaWUc7QIKKVUFNMioJRSUUyLgFJKRTEtAkopFcW0CCilVBTTIhDlRKRVRMY4neNcRKRcRIyIuO3Hr4jInQFq+xIR2evxuEpErgxE23Z7FSJyWaDa83OaIiKPikiziKwP5bSHw/u9UMGnRSCIRGSF/SGMdzqLL8aYFGPMwUC3KyJ3iciqQLfbxxhznTHmMT9yGBEZd4623jLGTAxELhFZJCIPeLU/1RizIhDtD8IC4Cqg2BhznucIEblQRNpEJNX7RSKyRUS+4FF0W71utwUypPf7E8j3op9pvee9UVoEgkZEyoFLAAPcFMTpuIPVdjQYwfOvDKgyxrR5jzDGrAVqgA96DheRacAUYLHH4Ax7RaHv9nQwQysHGGP0FoQb8P+A1cAvgH96jVsE/B5YCpwG3gTKPMYb4EvAQaAR+G/AZY+7y273l8AJ4AGsYv49oBo4DjwOpNvPv81uJ81+fB1QB+R6TGucR67fAq8ArfZ08oFfAc3AHmC2R85vAQfs/2EXcIs9fDLQDvTY7bTYw+OBnwGHgXp7HiT6mH8x9nMb7fz32Vnd9vgVwKft++PseXjSfv7T9vCV9mva7By3AZdhfQF+054PT/QN85h2FfBt+39qBh4FEjzm/yqvrMbOcA/QBXTa03vRo70rPebBr4Bj9u1XQLw9ri/b1+33sRa4e4BlrBB4wV4OKoHP2MM/5TX/7+/ntd8BlnsN+ynwnH2/3HN++7G8pwOP2JmPYi2XMUN9f7zei28A2+3nPQLkYS2jp4HXgUyP5//Nfl9P2u1PtYf7em8KgWeBBuAQ8CWPts4DNgKnsJbXXzj9vRKU7yqnA4zUm/2hvBeYay98eR7jFtkL8PvsL4UHPb9Y7A/GG0AWUArs450vvLuAbuCLgBtIBD5pT28MkAI8Bzzh0d6f7WlmY33xvN9rWp5FoNHOnAAstz8Yn8D6Un4AeMPjtR+yP0Qu+wPcBhR45PT+svwV1pdWFpAKvAj8p4/59zmsolNiP/8NfBeBxcB37RwJwIL+/j/78WX2/Psve94n0v8Xz06Paa8GHhjg//Kehw94ja/inSLwQ2AdMArIBdYAP/LK9kMgFrgeOIPHl5xXu29iFe0EYBbWF9kVvnJ6vbYEa7kstR+7sArQB+zH5QyuCPwd+AOQbP9v64HPDvH98X4v1mF98RdhFcfNwGz7/VsOfN/j+Z/EWrb6iu1Wr8/dAx6PXcAmrBW2OKzPz0HgGnv8WuDj9v0U4AKnv1eCcXM8wEi8YfXHdgE59uM9wFc9xi8C/uLxOAVrra3EfmyAaz3G3wsss+/fBRz2mt4y4F6PxxPt6fd9YWZgrX3vAP7g9VrvL7CHPcZ9Edjt8Xg69lq9j/97K3CzR07PwiZYRWKsx7ALgUM+2loOfM7j8dX4LgKPAw9h9X97t9Pfl0wn9pq9xzDvLx7PaV8PHOjv//IxDwcqAgeA6z3GXYPVbdOX4yweX7xYX3rv+fLB+hLvAVI9hv0nsMhXzn7aeB34jn3/KqwVgFj7cbn9f7V43Sb3004e0IHHVh1wO/YKwxDeH+/34mMej58Ffue1jP7dx/+XYbfft1X8rvcGOJ/3fpa+DTxq318J3I/9OR6pN90nEBx3Aq8ZYxrtx0/Zwzwd6btjjGnF2qQv7G88VjePr3HY46q9nu/G+nBijGnB2kyeBvz8HNnrPe6f7edxSt8DEfmEiGwVkRYRabHbz/HRbi6QBGzyeP6r9vD+FPLeeeDLv2MVmfX2kTifHOC5AA3GmPZzPGeg+T8c/b1Xnm03GWO6PR6fwWOee7Vzwhhz2qutokFkeQxrKw/g48BTxpgur+fkGGMyPG67+2mnDGvLpdbjvf0D1hYBDP798ebXMikiMSLyExE5ICKnsAoI+F4my4DCvsx27u9gf26wutUmAHtEZIOIvH+QuSPCSN0p5hgRSQQ+DMSISJ09OB7IEJGZxpht9rASj9ekYHU7HPNoqgSosO+Xeo0zXpM9hrVA9ynF6laot9ufhbWZvBj4H+DaIf1zHkSkDHgYuAJYa4zpEZGtWB/2/jI2Yn1gpxpjjvoxiVo85hHW/9QvY0wd8Bk71wLgdRFZaYyp9PUSP6bvPe2++d+GVcywp5c/yLb73itf762/jgFZIpLqUQhKsfrj/fUc8FsRWQj8G9Za+FAcwdoSyPEqYMCQ3p+h+ihwM3AlVgFIx9qn42uZPIK1JTq+v8aMMfuB20XEhTV/nhGRbNPPzvZIplsCgfcBrM30KVj9tLOwdpS+xTtrXQDXi8gCEYkDfgS8bYzxXPv8hohkikgJ8GVgoKMyFgNfFZHRdkH5D6ydb90ikgA8ibWGczdQJCL3BuD/TMb6UDUAiMjdWFsCfeqBYvv/wxjTi1U0fikio+zXFInINT7a/yvwJREpFpFMrJ3Q/RKRD4lIsf2w2c7V45FjKL+DuM+edhbWvOub/9uAqSIyy563P/B63bmmtxj4nojkikgOVn/0k4MNZy8ra4D/FJEEEZmBteb650G00QY8g7Xju9oYs3GwOex2aoHXgJ+LSJqIuERkrIhcCkF7f/qTilWMmrAK9X94jfee1nrglIh8U0QS7S2JaSIy3859h4jk2stui/2aHkYYLQKBdydWn+JhY0xd3w34DfAxj0MSnwK+j9UNNBf4mFc7/8DaabUVeAnrqAhf/oR1lMtKrB257Vh9pWD1E9cYY35njOkA7gAeEJF+1378ZYzZhdW1tBbrwzUdawdqn+VYa7t1ItLXLfZNrB3Y6+zN9dex9l/052FgCdaX7mastVZf5gNvi0gr1o7nLxtjDtnjfgA8Zm/uf3gQ/+JTWF9sB+3bAwDGmH1YO25fB/YD3r+FeASYYk/v7/20+wDWESfbsfbRbO5rewhux+q7PwY8j7WDdOkg23gMa8vkcR/jW7x+J/A1H8/7BNbO1b4jqp4BCuxxwXh/+vM4VpfYUTvHOq/x73pvjDE9wI1YK2qHsLZW/4i1BQHWFnOFnftB4CN+dCNGHLF3gKgQEpFFWF/M3/Mx3gDjg7C5rJRS76JbAkopFcW0CCilVBTT7iCllIpiuiWglFJRLCJ+J5CTk2PKy8udjqGUUhFl06ZNjcYYXz/IBCKkCJSXl7Nx45AOYVZKqaglIgP90h7Q7iCllIpqWgSUUiqKaRFQSqkopkVAKaWimBYBpZSKYloElFIqimkRUEqpKKZFQCmlolhE/FhMqZGmq6eXl3fUsqSijsbWTnJS4rh1bjGXTRiFyyXnbkCpANEtAaVCbFP1Ca791Ur++NZBijKSuHJyHoUZifz4pd3c9tBaGls7nI6ooohuCSgVQv/YepTvv1DBXReWc97oLET61vrTWThhFM9uruHGX6/i2c9fRGFGoqNZVXTQLQGlQuT5zTU88M9dfPu6yZw/JtujAFhcLuFD80q4YvIoPvGn9Zxq73IoqYomWgSUCoG3DzZx/z938Y1rJlGalTTgc6+fVsCYnGS+/vQ29HofKti0CCgVZCfaOrnvqc18/tKxlJyjAACICHdcUMa+46d5dlNNCBKqaKZFQKkgMsbwned3cMGYbGYUZ/j9utgYF5+/dCwPvLSb+lPtQUyoop0WAaWC6LVd9VQcPcmH5pYM+rVl2clcOjGX/3h5dxCSKWXRIqBUkHR29/LAS7u444Iy4txD+6jdPLOIVfsb2Xy4OcDplLJoEVAqSJ5cV0VuSvyguoG8JcbF8OF5Jfzon7t0J7EKCi0CSgXB2c4efvPGAW6bXzrsthaMy6HxdAcr9zcGIJlS76ZFQKkgeHrDYcaPSjnn4aD+cLmEW2YX8bMle3RrQAVc0IqAiJSIyBsisltEKkTky/bwLBFZKiL77b+ZwcqglBO6enr5w8qDvH9GQcDaPH9MNifPdrOqUrcGVGAFc0ugG/i6MWYycAFwn4hMAb4FLDPGjAeW2Y+VGjFe3lFLTko840alBqxNlwjXTcvnD28eCFibSkEQi4AxptYYs9m+fxrYDRQBNwOP2U97DPhAsDIo5YRFq6u4anJewNu9eFwOu2tPs6fuVMDbVtErJPsERKQcmA28DeQZY2rBKhTAKB+vuUdENorIxoaGhlDEVGrYdh07RU3zGeaUBb6XMzbGxZWT83h45cGAt62iV9CLgIikAM8CXzHG+L0KY4x5yBgzzxgzLzc3N3gBlQqgx9dWsXDSKGKCdE2AhZNGsaSinpNn9eRyKjCCWgREJBarAPzZGPOcPbheRArs8QXA8WBmUCpUznR289L2Wi6d0O/GbUCkJ8YysySd5zbrOYVUYATz6CABHgF2G2N+4THqBeBO+/6dwD+ClUGpUHp5Rx2TClLJSo4L6nQWThzFk+uq9XBRFRDB3BK4GPg4cLmIbLVv1wM/Aa4Skf3AVfZjpSLeX9YfZsG44HddTilIo72rl82HW4I+LTXyBe3KYsaYVYCvjtErgjVdpZxQ3dRG5fFWvnzF+KBPS0RYMD6HZzYdYW4QdkCr6KK/GFYqAJ7ZWMOFY7Nxx4TmI7VgXA4vba+lo7snJNNTI5cWAaWGqbfX8OzmGi4ZH7qj2HJS4inPTmb5bj2uQg2PFgGlhml91Qni3C7Ks4d/nqDBuGhcNn/TK4+pYdIioNQwPbOxhovH5bznwvHBNr88i7cPNnHyjP5mQA2dFgGlhuFMZzdLKuq4aGxOyKedFOdmenE6S3bVhXzaauTQIqDUMLy6s44J+cH/bYAv54/O5vktRx2ZthoZtAgoNQx/2XCES8aFfiugz+zSDHbUtNDY2uFYBhXZtAgoNURHTpxhb93poJwszl/x7hhml2by6k7tElJDo0VAqSH664YjXDQ2m9gQ/TbAl/llWbywVbuE1NBoEVBqCLp6evnLhiMsnBi8k8X5a2ZJBrtqT2mXkBoSLQJKDcGy3cfJTY2nJADXEB6uOLdLu4TUkGkRUGoInlhXxWUTw+c6F/PLsnhx2zGnY6gIpEVAqUE62NBKxdFTnD862+ko/zKzJIOdR0/SpF1CapC0CCg1SI+utq4eFucOn49PnNvFzJIMlu6qdzqKijDhsxQrFQFOnu3i71uPcmUQLiQ/XPPKsnhxu3YJqcHRIqDUIPx1w2FmFmc49gvhgcwuzWDL4RY9l5AaFC0CSvmpq6eXR1ZVce20fKej9CshNoZpRem8vlu7hJT/tAgo5aeXd9SSkxLH2NwUp6P4NLc0k5d31jodQ0UQLQJK+cEYw+/fPMB10wqcjjKgOaWZrD3QxJnObqejqAihRUApP7x96ASn27uZVZrhdJQBpSS4mZCXysp9DU5HURFCi4BSfnh45UGunpKPK8QXjhmK2aUZvLRDu4SUf7QIKHUO1U1tbKg+wSXjnTtl9GDMK8tixd4Gunp6nY6iIoAWAaXO4dHVVVw2YRQJsTFOR/FLVnIc+WkJbDh0wukoKgJoEVBqAO1dPTy/5ShXTHL+bKGDMas0g9f0spPKD1oElBrAS9trGZObzKi0BKejDMrc0kyW7jqOMcbpKCrMaRFQagBPrKtm4YTI2goAKM1Koru3l331rU5HUWFOi4BSPlQeb+XwiTPMLgvvw0L7IyLMLs3k9d3aJaQGpkVAKR+e31zDRWOzcbsi82MyszidZbuPOx1DhbnIXLqVCrLeXsPzW45y0djIOCy0P1MK0tlde5qTZ/WEcso3LQJK9WPT4WbcMS7Ks52/fORQxbldTClMY9X+RqejqDCmRUCpfjy3uYYLx2YjEfAL4YFML0pnmZ5VVA1Ai4BSXnp6DUsq6rkgjC4fOVSzSjJYsa9BDxVVPmkRUMrLxqoTZCTGkp8eWb8N6E9eWgLxbhd76k47HUWFKS0CSnl5aUct88oznY4RMNOL0nW/gPJJi4BSHnp7Da/urOO88sjvCuoztTCdN/bqoaKqf1oElPKwraaFeLeLosxEp6MEzNTCNLYcaaG9q8fpKCoMBa0IiMifROS4iOz0GPYDETkqIlvt2/XBmr5SQ7F0Vz1zy0ZOVxBAcrybsqwkNlY1Ox1FhaFgbgksAq7tZ/gvjTGz7NvLQZy+UoO2pKKOOaUjqwiAtTWwqlKvNqbeK2hFwBizEtATmquIUd3Uxom2TsaOCt8LyQ/VlII03tKdw6ofTuwT+IKIbLe7i0beKpeKWK9V1DG3LDMiLiE5WOPzUjnY0KankFDvEeoi8DtgLDALqAV+7uuJInKPiGwUkY0NDboZq4JvSUU9s0pG5npJbIyLifmprNerjSkvIS0Cxph6Y0yPMaYXeBg4b4DnPmSMmWeMmZebmxu6kCoqnTzTxa7aU0wrSnM6StBMLkhl1X5doVLvFtIiICIFHg9vAXb6eq5SobRi33GmFqYR746M6wgPxdTCdFZV6n4B9W7uYDUsIouBy4AcEakBvg9cJiKzAANUAZ8N1vSVGowlO+tGbFdQn9HZydSf6qDhdAe5qfFOx1FhImhFwBhzez+DHwnW9JQaqs7uXt6qbOS/PjjD6ShB5XIJUwvTWHuwiZtmFjodR4UJ/cWwinobqk5QkJ5AZlKc01GCblJ+qp5HSL2LFgEV9V6rqGP2CO8K6jO1MJ01B7QIqHdoEVBRzRjDa7vqmTPCThXhS3FmIq3t3dQ0n3E6igoTWgRUVNtTdxpjoGQEnTBuICLWfoE1B5qcjqLChBYBFdWW7qpnTllGxF9GcjAmF6Tx1j79vYCyaBFQUW1JFO0P6DO9KJ3VB5r0kpMK0CKgotixlrMcOXGGyQUj91fC/RllX3Jyb71eclJpEVBRbOmueuaUZhLjip6uoD7TCtP0UFEFaBFQUezlHbVRc1SQt6mF6azU/QIKLQIqSrWc6WTH0ZPMKE53OoojphSmsbG6mY5uveRktNMioKLSst3HmV6UPqJPGDeQ1IRYSjL1kpNKi4CKUi/tqB1x1xIerBnF6Szfc9zpGMphWgRU1Gnr6GbdwSZmj8BrCQ/GzJIMLQJKi4CKPm/ua2BiXiop8UE7iW5EGJ2TTHNbJ0dO6CkkopkWARV1XtoevUcFeXKJMKskgxV7dWsgmmkRUFGlo7uHN/c1ME+LAGB1Cb1aUed0DOUgLQIqqqw50ERJViIZUXDtAH/MKslgy+EWTp7tcjqKcohfRUBEnhWRG0REi4aKaK/uqGNOlO8Q9pQQG8PUwjTe0B3EUcvfL/XfAR8F9ovIT0RkUhAzKRUUPb2GpbvqmF+e5XSUsDK7NJOXd9Q6HUM5xK8iYIx53RjzMWAO1gXil4rIGhG5W0RigxlQqUDZcriZtKRY8tISnI4SVuaWZbL6QCNnO/XXw9HI7+4dEckG7gI+DWwBHsQqCkuDkkypAFtSUc9c7Qp6j7SEWMaPSmXZnnqnoygH+LtP4DngLSAJuNEYc5Mx5mljzBeBlGAGVCpQXtul+wN8uWBMNs9tPup0DOUAf38t80djzMueA0Qk3hjTYYyZF4RcSgXUgYZWWtu7GZ2T7HSUsDS/PJMn1lbRcqZTj5yKMv52Bz3Qz7C1gQyiVDAt3VXP3LLMqLqM5GAkxbmZWZLByzv0NwPRZsAiICL5IjIXSBSR2SIyx75dhtU1pFREWFJRx+zSDKdjhLULx2bz141HnI6hQuxc3UHXYO0MLgZ+4TH8NPCdIGVSKqBaznSyt+40X7ligtNRwtrskkweXV1F5fHTjBuV6nQcFSIDFgFjzGPAYyLyQWPMsyHKpFRArdzfyNTCNOLc+lvHgcS4hAXjcnh6Qw3fvWGy03FUiJyrO+gO+265iHzN+xaCfEoN27Ld9Uwv0q4gf1w6IZdnN9fQ1dPrdBQVIudaNeo7lCIFSO3nplRY6+k1rNzXoPsD/FSYkUhBegLLdutpJKLFubqD/mD/vT80cZQKrO01LaQnxpKTEu90lIhx6YRcnlxXzbXT8p2OokLA3x+L/VRE0kQkVkSWiUijR1eRUmFrxd4GphdF58Xkh+r80dlsP9qiF5uJEv7uKbvaGHMKeD9QA0wAvhG0VEoFyIq9x5lerF1BgxHndrFgXA6L1x92OooKAX+LQN9J4q4HFhtjTgQpj1IBc/JsF/uPtzIxT3dfDdZlE0bxt401dOsO4hHP3yLwoojsAeYBy0QkF2gPXiylhm9NZSOT8lP10NAhKMlKIjsljjf2NjgdRQWZv6eS/hZwITDPGNMFtAE3BzOYUsO1Ym8DUwt1f8BQvW9CLk+9Xe10DBVkg1lFmgzcJiKfAG4Frg5OJKWGzxjDW/t1p/BwXDgmm43VzdSd1I3+kczfo4OeAH4GLADm2zc9e6gKW4dPnKGju5fizESno0SshNgYzivP4vktNU5HUUHk76mk5wFTjDHG34ZF5E9YRxMdN8ZMs4dlAU8D5VhXKPuwMaZ5MIGV8seqykamFaXrWUOHacG4HJ5YV83nLh2r83KE8rc7aCcw2F+OLAKu9Rr2LWCZMWY8sMx+rFTAvbWvkSkFaU7HiHgT81M529XDjqMnnY6igsTfIpAD7BKRJSLyQt9toBcYY1YC3oeS3gw8Zt9/DPjAoNIq5YfeXsPag01M0/0BwyYiXDwuh79t1C6hkcrf7qAfBGh6ecaYWgBjTK2IjPL1RBG5B7gHoLS0NECTV9FgV+0p0hLdZCXrFbICYcG4HH74YgXfv3EK7hg93Hak8fcQ0Tex+vBj7fsbgM1BzIUx5iFjzDxjzLzc3NxgTkqNMKsrG5mqXUEBk5eWQHZKPOsO6m9ERyJ/jw76DPAM8Ad7UBHw9yFMr15ECuw2CwA9VaEKuFWVjUwu0K6gQDpvdBZ/36oXoh+J/N22uw+4GDgFYIzZD/jsyhnAC8Cd9v07gX8MoQ2lfOrq6WVTdTOTC/RUEYF0wZhsXquoo7NbTyMx0vhbBDqMMZ19D0TEDQx4uKiILMa6GP1EEakRkU8BPwGuEpH9wFX2Y6UCZntNCwXpCaQmxJ77ycpvOSnxFGUmsrqy0ekoKsD83TH8poh8B+uC81cB9wIvDvQCY8ztPkZdMYh8Sg3KmgNNTNb9AUExtyyTl3fUsnDSUDoBVLjyd0vgW0ADsAP4LPAy8L1ghVJqqFbtb2RyvhaBYJhXlsWy3fX09Pr9m1EVAfw9OqgXa0fwvcaYW40xDw/m18NKhUJ7Vw/bj55kku4PCIq8tAQykuLYWKVHCY0k57rQvIjID0SkEdgD7BWRBhH5f6GJp5T/thxuoSQziaQ4f3s51WDNKcvklZ11TsdQAXSuLYGvYB0VNN8Yk22MyQLOBy4Wka8GPZ1Sg7DmQCOT83UrIJjmlWWydFc92hEwcpyrCHwCuN0Yc6hvgDHmIHCHPU6psLG6slF3CgdZaVYSXT29VB5vdTqKCpBzFYFYY8x7jgkzxjTwziUnlXLc2c4edteeZqJuCQSViDC7NIPXd9c7HUUFyLmKQOcQxykVUpuqmynPSSIhNsbpKCPerJJMXtulRWCkOFcRmCkip/q5nQamhyKgUv5Yc6CRSXpoaEhMKUhjX91pmtt0PXAkGLAIGGNijDFp/dxSjTHaHaTCxqr9ev2AUIlzu5halM6KfXrqr5FAzwurIt7p9i72H29lQp7uDwiVGUXpLNutRWAk0CKgIt6GqhOMz0shzq2Lc6jMKsngrf2N+uvhEUA/NSrira5sYpIeFRRS2SnxZCbFsq2mxekoapi0CKiIt2p/I1ML9foBoTajOIM39miXUKTTIqAi2om2To40n2FMbrLTUaLOzOJ0lmsRiHhaBFREW3ugiSkFabhduiiH2oS8VKoa22hq7XA6ihoG/eSoiLZyX4OeKsIh7hgX04rSWaUXmoloWgRURFtV2cj0It0f4JSphdolFOm0CKiIdbjpDGc6uynOTHQ6StSaWZzOW/sa6NVDRSOWFgEVsVYfsLYCRMTpKFFrVFoCSfFuKo6dcjqKGiItAipirdh7nCmFuj/AaTOK9RQSkUyLgIpI3T29rDnQxIziDKejRL0ZRRm6XyCCaRFQEWnrkRZyU+PJTIpzOkrUm1yQxt6605w82+V0FDUEWgRURFqxt4EZelRQWIhzu5hckMaq/XqoaCTSIqAi0vK9x/XQ0DAyrTBdTyERobQIqIjT2NpBdWObnjo6jMwssXYO6wXoI48WARVxlu85zoySDNwxuviGi4L0ROLdMXqoaATST5GKOK9V1DNLjwoKOzNL9NfDkUiLgIoo7V09rD3YyKxSLQLhZmZxBq/v1gvQRxotAiqirD3YRFlWMmkJeonrcDO5II0Dx1v1rKIRRouAiihLdtYxq0S3AsJRrH1W0Tf3NTgdRQ2CFgEVMXp6Da9V1HHe6CynoygfZpVksKSizukYahC0CKiI8fahJjKT48hLS3A6ivJhTmkmqyobae/qcTqK8pMWARUx/rmtlvnluhUQztISYynLSmbtgSanoyg/aRFQEaGn1/BqRR3nj852Ooo6hzllGbyqXUIRQ4uAighrDzSRmRRLfrp2BYW7uaVZLN1VT49eaCYiaBFQEeFvm45w8bgcp2MoP+SnJ5CRGMvGqhNOR1F+0CKgwl5bRzfLdh/norFaBCLF/NFZvLjtmNMxlB8cKQIiUiUiO0Rkq4hsdCKDihyv7KxjckEq6Yn6A7FIcf7oLF7ZWaddQhHAyS2BhcaYWcaYeQ5mUBFg8frD2hUUYQrSE0lPimWDdgmFPe0OUmFtb91pqhrbmFuW6XQUNUjnj87i+S1HnY6hzsGpImCA10Rkk4jc098TROQeEdkoIhsbGvRn6NHq8bVVLJyYi9ul6yuR5uKxObyyo1Z/OBbmnPpkXWyMmQNcB9wnIu/zfoIx5iFjzDxjzLzc3JAG+pUAAA9PSURBVNzQJ1SOa+3o5oVtx1g4Kc/pKGoIslPiGZObwrLdenrpcOZIETDGHLP/HgeeB85zIocKb4vfrmZ6UTpZyXox+Uh10dhs/rrxiNMx1ABCXgREJFlEUvvuA1cDO0OdQ4W3ju4eHn7rEO+fUeh0FDUM88uz2Hy4mfpT7U5HUT44sSWQB6wSkW3AeuAlY8yrDuRQYezvW45SmJHI6Jxkp6OoYUiIjeHCMdksfvuw01GUDyEvAsaYg8aYmfZtqjHmx6HOoMJbZ3cvv15eyY0zdStgJLh80iieWn+Y7p5ep6OofughFyrs/PntakalxjOlIM3pKCoAyrKTyU6JY5lefzgsaRFQYaW1o5tfL6/kw/NKnI6iAuiKSXk8suqQ0zFUP7QIqLDyi9f2MrM4nbJs3Rcwkpw/JotDDa3sPHrS6SjKixYBFTYqjp3k+S1H+cj8UqejqABzu1xcPTWf3795wOkoyosWARUWunp6+fdntvOheSWk6YniRqTLJ41i5b4Gjpw443QU5UGLgAoL/7NsP/FuF5dN0F+Hj1RJcW4un5TH/75R6XQU5UGLgHLcxqoTPLmumk9fMgYRcTqOCqLrpuXz0o5aak+edTqKsmkRUI5qbO3gvj9v5tMLxpCZpKeHGOnSEmO5dEIuv31D9w2ECy0CyjE9vYYvPrWFC8fmMEdPFR013j+jkH9sPcrRFt0aCAdaBJRjfrZkL60d3dw6t9jpKCqE0hNjuXxSHg++vt/pKAotAsohr1XU8ezmGr6wcBwxLt0PEG1umFHAkoo6DjW2OR0l6mkRUCF3uOkM33x2O19YOE4PB41SKfFurpuWz09f3eN0lKinRUCFVEd3D597chM3zixkfF6q03GUg66Zms/6QyfYXtPidJSopkVAhdR/vbKHlAQ3107NdzqKclhCbAy3zC7iR//chTHG6ThRS4uACpkVe4/z4rZjfHrBaP09gALgsomjqDvVznI9w6hjtAiokGhq7eAbf9vOZy8dS2qC7gdQlhiX8JH5pTzw0m693oBDtAiooDPG8M1nt3Ph2GymFqY7HUeFmdklGaQmuHlyXbXTUaKSFgEVdIvXH+ZAQ5v+HkD1S0S44/wyHly2nxNtnU7HiTpaBFRQHWxo5aev7uXey8YSG6OLm+pfSVYSF4zJ1kNGHaCfShU07V093Pvnzdwyp4jizCSn46gw98E5xSzdVc/mw81OR4kqWgRU0Nz/YgWZyXFcNTnP6SgqAiTHu7n9vFK+9ex2unQncchoEVBB8dTb1azc16iHg6pBuWhsNinxbn6r1xwIGS0CKuDWVDby30v28vWrJpAU53Y6joogIsInLx7No6ur9HrEIaJFQAXU5sPN3Pvnzdy3cBwFGYlOx1ERKDslno+eX8qXFm/hTGe303FGPC0CKmA2VTfzqUUbuOd9Y/T3AGpYFozLoTQ7iW8/t0NPKRFkWgRUQCzbXc8n7QIwu1QvEKOGR0S466Jyth5pYdGaKqfjjGhaBNSwGGP43YpK/v2Z7Xz9qgnMKtECoAIj3h3D166cwK+XV7J8T73TcUYsLQJqyE61d/G5Jzfx3Oaj3H/TVD01tAq4UWkJfPmK8Xz16W1sqDrhdJwRSYuAGpLtNS3c8OBbGAPfu2EK2SnxTkdSI9SEvFTuvWwsn3l8I5uq9YdkgaZFQA1Kb6/h9ysO8IlH1nPL7GLuvng0cW5djFRwzSjO4LPvG8MnF23gzX0NTscZUfTTq/xW3dTGRx5aywvbj/HDm6dy4dhspyOpKDKrJJOvXjmBr/xlC4tWH9KjhgJEi4A6p47uHn77RiU3/WY1E/JT+e51k8lNTXA6lopCE/NT+f6NU1m0por7ntpMyxk96+hwaRFQPvX0Gp7fUsMVP3+T5XuOc/9NU7lheiEul54GQjknLy2B+2+aBsBVv1jJs5tq6O3VrYKhkkjYpJo3b57ZuHGj0zGiRv2pdp7bXMPja6vJTIrjltlFTCvSH3+p8LO//jRPvl2NMXDfwnFcNz2feHeM07HChohsMsbMG/A5WgSUMYZDjW28sbeBV3bUsrf+NPPLM7l8Uh5jc1OcjqfUgIwxbKs5yas7a6luOsN10/O5YXoh543OivqDFrQIqH4ZYzjacpYNVSdYXdnEmspGOrp7mVWSwazSDGYUZUT9h0dFpuOn2ll3qIlN1c0ca2lnwbgcrpqSx8JJo8hKjnM6XshpEVAANLd1sqv2FDuOnmRzdTNbj7TQ2dPL5II0JoxKZVpRGkUZiXrKZzWitJzpZMuRFrbXtLCj5iQT81O5blo+10wtoDQ7Oi5yFLZFQESuBR4EYoA/GmN+MtDztQj470RbJ9tqWth+pIVtNSepOHaS0+3djM5Jpiw7idE5KYzLTSEvLV6/9FXU6OzuZecxayVo0+FmcpLjuXzyKC4Zn8PcsswRe8rzsCwCIhID7AOuAmqADcDtxphdvl6jRaB/rR3d7K49xc6jJ9lkr+E3t3UydlQK5dlJlGenUJ6TRF5aAi79wlcKsH7wWNnQyvaaFvbUnabyeCslWUlMzk9lYn4qJVlJFGcmkpeWQE5KPAmxkbuj2Z8i4ET5Ow+oNMYcBBCRvwA3Az6LwHCcPNtF/al2jAGDsf4a6DWGXmPo6bX+dvcYunsNXT29/7rf3dtLT2/fc6y+9L6aaegbBj3G0NPTS1ePobOnl/auHtq7ejjb1UNbh3W/vbuHzu532u4rvi4R3C4h1u3CHSPEu2OId7uIi3ERGyO4XIIx0NHdS1tHN81tndSdaudoy1m6egz5aQn2Gn4yd19UTmFG4nu+8Fvb9ZzsSnnKT0sgf0o+V0/Jp6unlyPNZzhy4iybqpt5taKOptZOGls76OqxPqeJsTGkJLhJioshwR2DO0aIcQkC/9qiFsH+3LpIjIshOc56TUp8LCnxMaTEu0mKf6eNhNgYYmPsz75LcIl9c4HVstXmmJxk3DHB20fnRBEoAo54PK4Bzvd+kojcA9wDUFpaOuSJff7JTaw50DTk14eztAQ37V097K07zd6607zqdCClRpDkODcGaOvo5qy9UueEn946gw/PKwla+04Ugf76Jd7TJ2WMeQh4CKzuoKFO7KnPXDDUlyql1IjnxHGANYBnWSsGjjmQQymlop4TRWADMF5ERotIHPAR4AUHciilVNQLeXeQMaZbRL4ALME6RPRPxpiKUOdQSinlzD4BjDEvAy87MW2llFLv0HMDKKVUFNMioJRSUUyLgFJKRTEtAkopFcUi4iyiItIAVPv59BygMYhxgiHSMkdaXtDMoRBpeSHyMg82b5kxJnegJ0REERgMEdl4rhMmhZtIyxxpeUEzh0Kk5YXIyxyMvNodpJRSUUyLgFJKRbGRWAQecjrAEERa5kjLC5o5FCItL0Re5oDnHXH7BJRSSvlvJG4JKKWU8pMWAaWUimJhXwRE5FoR2SsilSLyrX7Gl4rIGyKyRUS2i8j19vByETkrIlvt2+89XjNXRHbYbf6PBPCK68PI+zGPrFtFpFdEZtnjVtht9o0bFai8fmYuE5Fldt4VIlLsMe5OEdlv3+70GO7kPO43r4jMEpG1IlJhj7vN4zWLROSQxzyeFai8w8lsj+vxyPWCx/DRIvK2Pe+ftk/N7mheEVnotRy3i8gH7HHBnsd/EpHjIrLTx3ixl8VKO/ccj3FOLMdDyhvw5di6bm543rBONX0AGAPEAduAKV7PeQj4vH1/ClBl3y8Hdvpodz1wIdZVzl4BrnM6r9dzpgMHPR6vAOY5OI//Btxp378ceMK+nwUctP9m2vczw2Ae+8o7ARhv3y8EaoEM+/Ei4NZwm8f241Yf7f4V+Ih9//d9y5XTeT2ekwWcAJKCPY/t9t8HzBngc3+9vSwKcAHwtlPL8TDzBnQ5DvctgX9dlN4Y0wn0XZTekwHS7PvpnOMqZSJSAKQZY9Yaa649DnwgzPLeDiwOUKZz8SfzFGCZff8Nj/HXAEuNMSeMMc3AUuDaMJjH/eY1xuwzxuy37x8DjgMD/prS6cy+2GuklwPP2IMeIwzmsZdbgVeMMWcClGtAxpiVWEXHl5uBx41lHZBhL6tOLMdDzhvo5Tjci0B/F6Uv8nrOD4A7RKQG6xoFX/QYN1qsbpc3ReQSjzZrztGmU3n73MZ7i8Cj9ubd/w3kJin+Zd4GfNC+fwuQKiLZA7zW6XnsK++/iMh5WGu5BzwG/9jevP6liMQHKG8gMieIyEYRWdfXtQJkAy3GmO4B2nQqb5+P8N7lOFjz2B8DLa+hXo79cc73IRDLcbgXAX8uSn87sMgYU4y1+fSEiLiwNpFKjTGzga8BT4lImp9tOpHXakDkfOCMMcazn/BjxpjpwCX27eMByutv5v8DXCoiW4BLgaNA9wCvdXoe+8prNWCt4T0B3G2M6bUHfxuYBMzH6hb4ZoDyBiJzqbFOFfBR4FciMtbPNp3K2zePp2NdQbBPMOexPwa7vAZzHvtjwOkHajkO9yLgz0XpP4XVN4oxZi2QAOQYYzqMMU328E1YlXKC3Waxx+sDeaH7Ief1GP+etSdjzFH772ngKazN9UA5Z2ZjzDFjzL/ZBfW79rCTA7zW0Xk8QF7sFYGXgO/Zm9h9r6m1N7s7gEcJn3nct8mPMeYg1v6h2VgnEcsQEbevNp3Ka/sw8LwxpsvjNcGcx/4YaHkN9XLsD5/vQ0CX48HsQAj1DevylweB0byzg2qq13NeAe6y70+2Z5Jg9ZHF2MPHYK2pZNmPN2DtaOnb2XO903ntxy77jR/j1WaOfT8Wqw/4cyGexzmAy77/Y+CH9v0s4BDWzrRM+344zGNfeeOw+rG/0k+7BfZfAX4F/CRM5nEmEO/xnP3YO2mxds567hi+1+m8HuPXAQtDNY89plGO7x2tN/DuHa3rnVqOh5k3oMtxQN+AYNywukz2Ya3Jf9ce9kPgJvv+FGC1vaBuBa62h38QqLCHbwZu9GhzHrDTbvM32F/CTua1x10GrPNqLxnYBGy3/58HsYtbCDPfivXlsw/4I/aXkj3uk0Clfbs7TOZxv3mBO4Aue7733WbZ45YDO+zMTwIp4TCPgYvsXNvsv5/yaHMM1tErlVgFId7pvPa4cqyVLpdXm8Gex4uxuoG7sFamPgV8DnulCeuL8X/t/2kHHkfcObQcDylvoJdjPW2EUkpFsXDfJ6CUUiqItAgopVQU0yKglFJRTIuAUkpFMS0CSikVxbQIKKVUFNMioJRSUez/A/d7kc83ZhDCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we use approximate bootstrap to calculate confidence\n",
    "# intervals for vaccine efficacy via Monte Carlo draws\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV  + np.random.normal(0, 1, B) * np.sqrt(var_rv)\n",
    "RUs = RU  + np.random.normal(0, 1, B) * np.sqrt(var_ru)\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sgE4OHh5zhy-",
    "papermill": {
     "duration": 0.009208,
     "end_time": "2021-01-04T05:11:00.895842",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.886634",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the code cell below  we calculate the effectiveness of the vaccine for the two groups that are 65 or older"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "pvzbZxlAzhy_",
    "outputId": "331589ef-abf0-476b-81af-40540e9c460f",
    "papermill": {
     "duration": 0.477817,
     "end_time": "2021-01-04T05:11:01.383093",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.905276",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall vaccine efficacy is 0.9471\n",
      "Two-sided 95% confidence interval of VE is [0.8134, 1.0518]\n",
      "One-sided 95% confidence interval of VE is [0.8438, 1]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEICAYAAABVv+9nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxjZ3no8d8j7/K+zni8jGcyM5mZ7MkACRAIexIotLeXQsoWSgktS4Hby3r5FNrmUyi9QOgGhAAJSwJtCCFwAyErCdk9+74v9tgeezzePd6k5/5xjhJFsWzZ1tHR8nw/H31s6Ryd99GR9OjV+756X1FVjDHGZJ+A3wEYY4zxhiV4Y4zJUpbgjTEmS1mCN8aYLGUJ3hhjspQleGOMyVKW4HOEiIyKyGq/45iPiLSJiIpIvnv9NyLyviQd+0oR2R91/ZiIvD4Zx3aPt1tErkrW8RIsU0TkByIyICLPpLLspYh9Low3LMEngYg84r7BivyOJR5VLVPVI8k+rohcLyJ/SPZxI1T1GlW9LYE4VETWzHOsx1T13GTEJSK3isiNMcc/T1UfScbxF+CVwBuAZlV9afQGEblCRMZEpDz2TiKyVUQ+GvWBOhpzeUcyg4x9fpL5XMxS1ouem1xlCX6JRKQNuBJQ4K0elpPv1bFzQRafv5XAMVUdi92gqk8CncCfRt8uIucDG4E7om6ucisBkcvPvAzapIiq2mUJF+DvgMeBrwO/jtl2K/Bt4H5gBPg9sDJquwJ/AxwBTgP/AgTcbde7x/0GcAa4EecD+QvAcaAX+CFQ6e7/Dvc4Fe71a4AeoD6qrDVRcf0n8Btg1C1nOXATMADsAy6JivOzwGH3MewB/sS9fQMwAYTc4wy6txcB/xc4AZxyz0FJnPOX5+572o3/I26s+e72R4C/dP9f457DIXf/n7m3P+reZ8yN4x3AVTjJ7TPuefhR5Laoso8Bn3Mf0wDwA6A46vz/ISZWdWO4AZgGptzyfhV1vNdHnYObgC73chNQ5G6LxPa37vPYDbx/jtfYCuAe93VwCPige/sHYs7/389y388DD8Xc9lXgLvf/tujzncDrvRL4nhvzSZzXZd5in5+Y5+JTwA53v+8By3BeoyPAA0B11P7/7T6vQ+7xz3Nvj/fcrAB+DvQBR4G/iTrWS4F2YBjn9fp1v/NK0vKT3wFk+sV9w30YuMx9YS2L2nar++J8lfuG/2Z00nBf9A8DNUArcIDnk9n1wAzwMSAfKAH+wi1vNVAG3AX8KOp4P3HLrMVJKm+JKSs6wZ92Yy4GHnJf9O/FSbg3Ag9H3fft7hsk4L45x4DGqDhjE+FNOAmpBigHfgV8Oc75+yucD5QWd/+HiZ/g7wD+jxtHMfDK2R6fe/0q9/z9s3vuS5g9qeyKKvtx4MY5HlfsObwxZvsxnk/w/wA8BTQA9cATwD/GxPYPQAFwLTBOVAKLOe7vcT6Qi4GLcZLU6+LFGXPfFpzXZat7PYDz4fLH7vU2Fpbg7wa+A5S6j+0Z4EOLfH5in4uncJJ6E84H3xbgEvf5ewj4YtT+f4Hz2op8kG6Led/dGHU9AGzGqYwV4rx/jgBvcrc/CbzH/b8MuNzvvJKsi+8BZPIFp/1zGqhzr+8DPhm1/Vbgp1HXy3BqWy3udQWujtr+YeBB9//rgRMx5T0IfDjq+rlu+ZFkWIVTa94JfCfmvrHJ6btR2z4G7I26fgFubTzO494GvC0qzugPLcH5ADgn6rYrgKNxjvUQ8FdR199I/AT/Q+BmnPbm2OPMlkCmcGvkUbfFJpXosq8FDs/2uOKcw7kS/GHg2qhtb8JpSonEcZaopIqT0F6UWHASdAgoj7rty8Ct8eKc5RgPAJ93/38Dzod7gXu9zX1cgzGXDbMcZxkwSdS3MeA63MrAIp6f2OfiXVHXfw58K+Y1enecx1flHj/ybfYFzw3wMl78Xvoc8AP3/0eBv8d9H2fTxdrgl+Z9wO9U9bR7/Xb3tmgdkX9UdRTna/aK2bbjNL3E24a77XjM/vk4bzxUdRDnq+v5wNfmif1U1P9nZ7leFrkiIu8VkW0iMigig+7x6+Ictx4IApuj9v+te/tsVvDicxDPp3E+QJ5xR6z8xRz7AvSp6sQ8+8x1/pditucq+tj9qjoTdX2cqHMec5wzqjoSc6ymBcRyG863M4D3ALer6nTMPnWqWhV12TvLcVbifOPojnpuv4NTk4eFPz+xEnpNikieiHxFRA6LyDDOhwPEf02uBFZEYnbj/jzu+wanqWsdsE9EnhWRtyww7rSVrR1PnhOREuDPgDwR6XFvLgKqROQiVd3u3tYSdZ8ynKaArqhDtQC73f9bY7ZpTLFdOC/WiFacr/qn3ONfjPPV9Q7gX4GrF/XgoojISuC7wOuAJ1U1JCLbcN7Is8V4GufNeJ6qnkygiG6izhHOY5qVqvYAH3TjeiXwgIg8qqqH4t0lgfJjy46c/zGcDyrc8pYv8NiR5yrec5uoLqBGRMqjknwrTvt3ou4C/lNEXgP8D5za82J04NTg62I+nIBFPT+L9efA24DX4yT3Spw+lHivyQ6cb5BrZzuYqh4ErhORAM75uVNEanWWjutMYzX4xftjnK/OG3HaRS/G6XR8jOdrSwDXisgrRaQQ+EfgaVWNrjV+SkSqRaQF+Dgw1+iFO4BPisgq98Pin3A6smZEpBj4MU7N5P1Ak4h8OAmPsxTnDdMHICLvx6nBR5wCmt3Hh6qGcT4QviEiDe59mkTkTXGO/1/A34hIs4hU43TozkpE3i4ize7VATeuUFQcixnn/xG37Bqccxc5/9uB80TkYvfcfinmfvOVdwfwBRGpF5E6nPbfHy80OPe18gTwZREpFpELcWqcP1nAMcaAO3E6kY+ravtC43CP0w38DviaiFSISEBEzhGRV4Nnz89synE+aPpxPoT/KWZ7bFnPAMMi8hkRKXG/AZwvIi9x4363iNS7r91B9z4hsoAl+MV7H04b3glV7YlcgH8H3hU1LO924Is4TTOXAe+KOc4vcTqAtgH/D2f0QDzfxxkN8ihOp+gETtskOO2ynar6LVWdBN4N3Cgis9ZaEqWqe3Cae57EeeNcgNMZGfEQTi21R0QiTVWfwekMfsr9Cv0ATn/BbL4L3IeTULfg1DbjeQnwtIiM4nTiflxVj7rbvgTc5n4F/7MFPMTbcZLWEfdyI4CqHsDpBH0AOAjEjvX/HrDRLe/uWY57I87IjB04fSJbIsdehOtw2sq7gF/gdDbev8Bj3IbzjeKHcbYPxoyD/19x9nsvTkdlZOTRnUCju82L52c2P8RppjrpxvFUzPYXPDeqGgL+CKcSdhTnW+YtODV/cL7p7nbj/ibwzgSa9jKCuJ0MxgMicitO0v1CnO0KrPXgK6wxxlgN3hhjspUleGOMyVLWRGOMMVnKavDGGJOl0mocfF1dnba1tfkdhjHGZIzNmzefVtVZf0iYVgm+ra2N9vZFDdE1xpicJCJxf/1tTTTGGJOlLMEbY0yWsgRvjDFZyhK8McZkKUvwxhiTpSzBG2NMlrIEb4wxWcoSvDHGZClL8MaYRTszNsXHf7qV/3j4IL3DWTGFelaxBG+MWZTekQne/u0nmJgO0X5sgHd+9ymmZsJ+h2WiWII3xizKl+/dx/rlFVz/8lV88MrV1AQL+fYjh/0Oy0SxBG+MWbCuwbM8sPcUf3xJEwAiwvte3sb3Hj9K1+BZn6MzEZbgjTEL9t3HjvDqdfWUFT0/X2FdWRFXrK7hZ892zHFPk0qeJXgROVdEtkVdhkXkE16VZ4xJjYnpEHe2d3L1ectftO3V5zbws2c7CIVtIaF04FmCV9X9qnqxql4MXAaM46wIb4zJYI8dPE1bXSm1ZUUv2tZWW0pZUT6PHzrtQ2QmVqqaaF4HHFbVuPMWG2Myw292dnNpa1Xc7Veuq+Onz55IYUQmnlQl+HcCd6SoLGOMR2ZCYR7c18tlK2vi7vPSthoePXCayZlQCiMzs/E8wYtIIfBW4L/jbL9BRNpFpL2vr8/rcIwxS9B+fIC6skLqy1/cPBNRFSykubqEp4+cSWFkZjapqMFfA2xR1VOzbVTVm1V1k6puqq+fdVlBY0yaeGDvKS5pid88E3FJaxX37e5JQURmLqlI8NdhzTPGZIUnDvezcUXlvPtd2lrN/XtOoWqjafzkaYIXkSDwBuAuL8sxxnhvZGKao31jnFNfNu++TVUl5OcJu7uGUxCZicfTBK+q46paq6pDXpZjjPFe+7EB1jSUUZg/f9oQES5oquQPNlzSV/ZLVmNMQp480s+6ZfPX3iPOW1HJowds4ISfLMEbYxLy5OF+NjZWJLz/xsYKtnYMMjFtwyX9YgneGDOv8akZDvaOsKahPOH7lBbl01odZMuJAQ8jM3OxBG+MmdeOziFW1pQm1P4ebeOKCv5w0Nrh/WIJ3hgzr20dg6yuL13w/c5bUcFjluB9YwneGDOvzccHEhoeGWttQzkHe0cYn5rxICozH0vwxph5be8YZE3DwhN8YX6AttpStp0Y9CAqMx9L8MaYOXUPnWVyJkzDHPPPzOXcZeU8c9TmpfGDJXhjzJy2nRhkbUMZIrKo+69bXs6TR/qTHJVJhCV4Y8yctpwYZNUiOlgjzl1Wzs6TQ0yHwkmMyiTCErwxZk7bOgZYXbfw9veI0qJ8llUU27w0PrAEb4yJKxxW9nQPs7pu8TV4gLUNZWw5bj94SjVL8MaYuI71j1FWlE9FScGSjrO6voz249bRmmqW4I0xce08OcTqRYx/j7W2oYxtHTapbKpZgjfGxLW9Y5C22uCSj7O8spiRiWn6RiaTEJVJlCV4Y0xc2zoGWbWEDtaIgIhbi7cfPKWSJXhjzKxCYWVfzwirltjBGrGqrtQ6WlPMErwxZlZH+kapLCmgrCg/Kcdb01BuHa0p5vWarFUicqeI7BORvSJyhZflGWOSZ0fnEOcs4QdOsVbXl7Kne5hw2BbiThWva/DfBH6rquuBi4C9HpdnjEmS7Z2DrKxNXoKvKC6gtDCfE2fGk3ZMMzfPEryIVACvAr4HoKpTqmo9LMZkiO0dg0kZIhltdX0pO0/acMlU8bIGvxroA34gIltF5BYRSV51wBjjmelQmP2nRpIyRDJaa02QHZ2W4FPFywSfD1wKfEtVLwHGgM/G7iQiN4hIu4i09/XZCuzGpIODp0apLy8iWJicDtaIVXWl7Oi0L/Kp4mWC7wQ6VfVp9/qdOAn/BVT1ZlXdpKqb6uvrPQzHGJOonScHlzTBWDxttU5Hq6p1tKaCZwleVXuADhE5173pdcAer8ozxiTPto5BVia5eQagKlhIUX6AjjNnk35s82Jej6L5GPATEdkBXAz8k8flGWOSYFvH4KLWYE3E6voy62hNkeQ2sMVQ1W3AJi/LMMYk18R0iCN9Y7QlcYhktKaqEvb1DPPmCxs9Ob55nv2S1RjzAnu6h2muLqEw35v00FoTtMU/UsQSvDHmBXZ42DwD0FITZH/PiGfHN8+zBG+MeYEtJwY9a54BWF5RTP/YJGOTM56VYRyW4I0xL7C9Y5BzGryrwecFhObqIPtPWS3ea5bgjTHPGTo7Te/IJM1VJZ6W01JdYs00KWAJ3hjznB2dg6yuLyUQEE/LaaoKssc6Wj1nCd4Y85ytJwZZnaQFPubSWhtkb7cleK9ZgjfGPGfz8QFPR9BEtFSXcLB31PNycp0leGMMAKrKjs5B1njYwRpRWVJAKKz0j9oi3F6yBG+MAaBz4CwBEWpKCz0vS0RoqS7hkNXiPWUJ3hgDwNYOp/Yu4m0Ha0RjVQmH+izBe8kSvDEGgK0nBlidxDVY59NYWcxBGwvvKUvwxhgAthwf8GQO+HhWVJVw4JTV4L1kCd4Y89wSfamswTdVlXDYmmg8ZQneGMO+7hGWVRQnfYm+udSXFTE4Pm1z0njIErwxhm2dg5yTwto7QCAgVov3mCV4Ywxbjp9Jaft7xIoqGyrpJUvwxhi2nvB2Bsl4llUUWQ3eQ5bgjclxwxPTnBqeoKU6+Ytsz6exsoTDvWMpLzdXeNqjIiLHgBEgBMyoqq3Pakya2dk5xOr6MvI8nkFyNo2VxfxuT0/Ky80Vqegyf42qnk5BOcaYRdjWMcCqFMwgOZvllcWcODNOOKyeT1Gci6yJxpgctyVFUwTPJliYT2lhPj3DE76Un+28TvAK/E5ENovIDbPtICI3iEi7iLT39fV5HI4xJlakicYvK6qKOXra2uG94HWCf4WqXgpcA3xERF4Vu4Oq3qyqm1R1U319vcfhGGOinRqeYGI6REN5kW8xLKso5ogleE94muBVtcv92wv8Anipl+UZYxZme8cgaxrKUzaD5GyWVRRzuNcmHfOCZwleREpFpDzyP/BGYJdX5RljFm5bxyBtdakfHhmtsbKEw31Wg/eCl6NolgG/cGsG+cDtqvpbD8szxizQ1hODvGJNna8xNFYWc8QSvCc8S/CqegS4yKvjG2OWRlXZ3TXEe65Y6WscDRVF9I1MMjUTpjDfBvYlk51NY3JU58BZCvMDVAe9X6JvLvmBAHXlhZw4M+5rHNnIErwxOWpH5xDn+Dg8MlpjZQnHbCRN0lmCNyZHbe8cZGWtvx2sEcsqimwsvAcswRuTo7Z1DPo2RUGsZRXFNm2wByzBG5ODVJU9XcOs8mEO+NksryjmyGlL8MlmCd6YHHS8f5xgYR6VJQV+hwI4bfDH+62TNdkswRuTg3aeHErpAtvzqS0rZPDsNONTtj5rMlmCNyYH7Tw5xMqa9EnwAREaK4s5dtpq8clkCd6YHLSj0/8pCmI1VtqskslmCd6YHBPpYG2rTZ8aPMCycutoTTZL8MbkmJODZynIC1Dl8y9YYy2vLObQKUvwyWQJ3pgcs+vkcNqMf4+2oqqEw32W4JPJErwxOWbXySFa0+QXrNEaK4s52j+GqvodStawBG9MjtneOZh27e8A5cUF5InQNzrpdyhZI6EELyI/F5E3i4h9IBiT4fZ2p18Ha0RTddDmhk+iRBP2t4A/Bw6KyFdEZL2HMRljPNI3MsnkTJi6svTqYI2wxT+SK6EEr6oPqOq7gEuBY8D9IvKEiLxfRNLjt87GmHnt7XY6WP1cg3UuzqRjtj5rsiTc5CIitcD1wF8CW4Fv4iT8+z2JzBiTdHu6h2mpTr8O1ogVlcUctFklkyahJftE5C5gPfAj4I9Utdvd9DMRaZ/nvnlAO3BSVd+ylGCNMUuzs3OI1pr0TfCNVSX2a9YkSnRN1ltU9d7oG0SkSFUnVXXTPPf9OLAXqFhMgMaY5NnTPcyVa/1dZHsuyyqK6B2eZGI6RHFBnt/hZLxEm2hunOW2J+e7k4g0A28GbllIUMaY5Ds7FaJr8CxNVSV+hxJXfiBgc9Ik0Zw1eBFZDjQBJSJyCRDpmakAEvmedxPwaaB8jjJuAG4AaG1tTeCQxpjF2H9qhObqEvLz0nu0c1N1CQd7R9nQaF/6l2q+Jpo34XSsNgNfj7p9BPj8XHcUkbcAvaq6WUSuirefqt4M3AywadMm+wmbMR7Z0zVMSxq3v0c0VhZz8JSNpEmGORO8qt4G3CYif6qqP1/gsV8BvFVErgWKgQoR+bGqvnuRsRpjlmBP11Baj6CJaKoqYW/PsN9hZIX5mmjerao/BtpE5H/FblfVr89yt8i2zwGfc49zFfC/Lbkb45893cNcc36j32HMq6k6yD3bu/wOIyvM1xgX+T1zGU47euzFGJMBVJUDp0bTeohkRGNlMV2DE0zNhP0OJePN10TzHffv3y+lEFV9BHhkKccwxizeycGzFOUHqEiTRbbnUpAXYFlFEcf6x1i3zOqRS5HoZGNfFZEKESkQkQdF5LSIWHOLMRliX/cIK9NwiuB4mqpLOGiLfyxZouOl3qiqw8BbgE5gHfApz6IyxiTV3p5hmqvTd/x7rMbKEg6cso7WpUo0wUe+110L3KGqZzyKxxjjgd1dw7TUpOcUwbNpri5hb7cNlVyqRBP8r0RkH7AJeFBE6oEJ78IyxiTTvu5hWjKoBt9SHWRfjyX4pUp0uuDPAlcAm1R1GhgD3uZlYMaY5JiYDtE1NJHWUxTEaqwqpmd4gonpkN+hZLREJxsD2IAzHj76Pj9McjzGmCQ71DvKisritJ+iIFp+IEBTldPRekFzpd/hZKxEpwv+EXAOsA2IfKQqluCNSXv7ekYyYoqCWM3VJezrGbYEvwSJ1uA3ARvVljs3JuPs7R5mRQY1z0Q0VZVYO/wSJfqdbRew3MtAjDHe2Ns9TGsGzEETq6UmyJ4uGyq5FInW4OuAPSLyDDAZuVFV3+pJVMaYpDlwaoR3viTzpuJuqQ5ywGaVXJJEE/yXvAzCGOONgbEpzk6FqCsr9DuUBasrK2RiOkT/6CS1ZUV+h5OREh0m+XvgGFDg/v8ssMXDuIwxSbCvZ4TW2lJEZP6d04yI0FZXyn5rh1+0ROei+SBwJ/Ad96Ym4G6vgjLGJMf+DJuiIFZzdQl7LcEvWqKdrB/BWcBjGEBVDwINXgVljEmOPd0jGZ7gg+zpGvI7jIyVaIKfVNWpyBX3x042ZNKYNJepI2giWmuCVoNfgkQT/O9F5PM4i2+/Afhv4FfehWWMWapwWDncN5qRP3KKaKkOcqRvlFDY6pOLkWiC/yzQB+wEPgTcC3zBq6CMMUvXMTBOaVE+pUULmZEkvZQU5lEdLORY/5jfoWSkhJ55VQ2LyN3A3ara53FMxpgk2NczwsoMrr1HtNYE2dc9wjn1ZX6HknHmrMGL40sichrYB+wXkT4R+bv5DiwixSLyjIhsF5HdIrKkZf+MMQuzr3uYpgzuYI1w5oa3X7QuxnxNNJ/AGT3zElWtVdUa4GXAK0Tkk/PcdxJ4rapeBFwMXC0ily85YmNMQvZ0D9OSwR2sES3VQfZYgl+U+RL8e4HrVPVo5AZVPQK8290WlzoiiyoWuBfrKTEmRfZ1Z+YskrFaamzxj8WaL8EXqOrp2Bvddvh5l2cXkTwR2Qb0Aver6tOz7HODiLSLSHtfnzXvG5MM41Mz9AxPsKKq2O9Qlmx5RTFnxiYZnZzxO5SMM1+Cn1rkNgBUNaSqFwPNwEtF5PxZ9rlZVTep6qb6+vr5DmmMScD+HucHTvmBzFnkI55AQGipDtqUBYsw37N/kYgMz3IZAS5ItBBVHQQeAa5eQqzGmATt7R6hNQuaZyKaq0sswS/CnMMkVTVvsQd2F+aeVtVBESkBXg/882KPZ4xJ3K6TQ1nR/h7RXB20kTSL4OX3t0bgYRHZgTP75P2q+msPyzPGuPZ0DbOyttTvMJKmpcZG0iyGZz9xU9UdwCVeHd8YM7twWDnQm11NNK01QQ6eGkFVM3LqY79kfg+MMeYFTpwZp7w4n7IMnqIgVmVJAXkBoWd4wu9QMooleGOyzJ7u7GqeiWitsZE0C2UJ3pgss+vkEC1ZMEVBrCYbKrlgluCNyTLbOwdZVZd9E3M1V5Wwp8s6WhfCErwxWURV2X1ymFV12ddEY1MWLJwleGOySNfQBIGAUFNa6HcoSddcXcKx/jFmQmG/Q8kYluCNySI7OwdZnYW1d4DigjxqS23xj4WwBG9MFtnROcTK2uwZ/x6rpSbI/p7R+Xc0gCV4Y7JKtnawRjRVl7CvxzpaE2UJ3pgskc0drBEt1UF220iahFmCNyZLdJw5S35ednawRrTaSJoFsQRvTJbY2jHA2mXlfofhqcjiHyMT036HkhEswRuTJbYcH2BVFk5REC0QEJuyYAEswRuTJTYfH2BtQ/Z2sEa01gTZawk+IZbgjckCkzMhDvWNsqo+u2vwAE1VQXafHPI7jIxgCd6YLLC7a5jm6iBF+YtehC1jrKy11Z0SZQnemCyw9cQg52Tx8MhoLTVBDvaOEg6r36GkPUvwxmSBZ472sybLR9BElBU5i5mcODPudyhpz7MELyItIvKwiOwVkd0i8nGvyjIml6kqzx4b4NwcSfAAq+pK2dVl7fDz8bIGPwP8rapuAC4HPiIiGz0sz5icdLx/nDwR6sqy9wdOsVbWBtnRaQl+Pp4leFXtVtUt7v8jwF6gyavyjMlVzx47w/rG8pxajHpVXSnbOwb9DiPtpaQNXkTagEuAp2fZdoOItItIe19fXyrCMSarPH3kDGtyYPx7tFV1ZeztHkbVOlrn4nmCF5Ey4OfAJ1T1RWObVPVmVd2kqpvq6+u9DseYrPPMsTM51f4OUFlSQFF+Hh1nzvodSlrzNMGLSAFOcv+Jqt7lZVnG5KLekQkGxqdoqcneOeDjWV1fyo6T1kwzFy9H0QjwPWCvqn7dq3KMyWVPHu5nY2MFgRxqf49otY7WeXlZg38F8B7gtSKyzb1c62F5xuScxw/1s355bjXPRKyuK2PriQG/w0hrXo6i+YOqiqpeqKoXu5d7vSrPmFz0xOHTbFxR6XcYvljbUMburmGmbRHuuOyXrMZkqJODZxmdmKG5usTvUHxRWpRPfXkR+7ptZsl4LMEbk6GePNzPxhW52f4esW5ZGe3Hz/gdRtqyBG9MhnrsQB/rl1f4HYavzqkv55mjluDjsQRvTAYKh5XHDvZxUXNutr9HnLusnM3HB+wHT3FYgjcmA+3pHiZYlE9DRbHfofhqWUURU6EwJwftB0+zsQRvTAZ6ZH8vFzTldu0dQETY0FjBU0esmWY2luCNyUAP7evlwhxvnonYsLyCxw7aPFazsQRvTIYZnphmb/cIGxpzu4M14oKmSh4/dNra4WdhCd6YDPPwvl7OW1GRE+uvJmJZRREBEQ73jfodStqxBG9Mhvntrh4uaa32O4y0ISKc31TB44f6/Q4l7ViCNyaDTM6EeOzgaS5trfI7lLSysbGSRw70+h1G2rEEb0wGeeJwP601QaqCubM8XyIuaKrkmaNnmJwJ+R1KWrEEb0wGuXdHN5dY7f1FKkoKaKkO2q9aY1iCNyZDTM6EuG93D5evrvU7lLR0YXMlD+495XcYacUSvDEZ4pH9fbTWBqkrK/I7lLR0cUs1D+61dvholuCNyRA/39xptfc5tNUGGZ8KccSGSz7HErwxGWDo7DSPHzrNy1ZZgo9HRLh0ZTX37e7xO5S0YQnemAxw15ZOLmqpoqwo3+9Q0tqmle+KP2AAAAtCSURBVNXcu9MSfISXi25/X0R6RWSXV2UYkwtUlR8+cZzXb1jmdyhpb+OKCo73j9E9ZLNLgrc1+FuBqz08vjE54ckj/YTRnF1ceyHyAwEuba3mvl1WiwdvF91+FLBBqcYs0ff/cJTXrm9AcnhpvoW4rK2aX+3o9juMtOB7G7yI3CAi7SLS3tdnU34aE21v9zBbTgzy6nX1foeSMS5squLAqRG6bBEQ/xO8qt6sqptUdVN9vb2IjYn2jfsP8OYLGm3myAUozA/wslU1/HLbSb9D8Z3vCd4YM7tdJ4doPz7A6zY0+B1KxrlidS2/2GoJ3hK8MWkoHFY+/4udvP2yZqu9L8L6xgoGxqbY1zPsdyi+8nKY5B3Ak8C5ItIpIh/wqixjss1Pnz3B9EyYV1nb+6IERHjl2nruePqE36H4ystRNNepaqOqFqhqs6p+z6uyjMkmR0+P8dXf7ud9L28jYCNnFu216xv4xdaTjE/N+B2Kb6yJxpg0MjEd4q9/vJk/ubSJlbWlfoeT0erKili/vIJ7tnX5HYpvLMEbkybCYeUTP9tGbVkRb7BfrSbFa9c38IPHj+XsgtyW4I1JA6rKF+/ZTeeZcW64crX9qClJLmiuZCoU5qF9uTmNsCV4Y3ymqnzlN/t44vBpPvH6dRTm29syWQIivO2iFXz9/gM5WYu3V5IxPvv3hw7x2909fPrq9ZTabJFJ95JVNYxOzuTkYiCW4I3x0U+eOs7tz5zgM1evp6K4wO9wslJAhHe+pJUv/Wo3E9O5tSi3JXhjfPLogT6+dv8BPv2m9VQHC/0OJ6td3FJFa02Qf3vokN+hpJQleGN8cKRvlI//dCsffc0allcW+x1OTnjXy1Zy+9PHaT+WO5PcWoI3JsVGJqb5wG3t/OllzWxorPA7nJxRU1rIX165mg//ZAu9IxN+h5MSluCNSaFQWPn4T7expqGM1623se6pdmlrNa9eV8/7f/Aso5PZ/wtXS/DGpIiq8oW7d9I3Msl7L1/pdzg5608uaWJFVQkfuPXZrO90tQRvTArMhMJ86Z7dPH3kDJ94/Vry8+yt5xcR4for2ijMD/CB27I7ydurzBiPHTg1wrtueZrtnUN87poNBAttrLvfAgHhQ686B1X40I82MzUT9jskT1iCN8YDp0cnuXNzJ+/7/jO84ztPsm5ZOZ9647mUFVtyTxd5AeGvrzqHs9MhPnr7FqZD2ZfkJZ1+vrtp0yZtb2/3OwxjFqVzYJxfbe/i3p09HDk9ygVNlVy2soZNK6spLrBFO9LVdCjMTQ8coKGimH+/7pKMaz4Tkc2qumnWbZbgjVm8yZkQ9+0+xY+fOs6+nmEuX1XLprYaNiwvz7hEkcumZsLc9OABaksL+bfrLqWkMHM+kC3BG5NEqsquk8P8fEsnv9x2ktaaIFed28BlK6spsKSesaZDYW75wxHOjE3xH39+acbMx28J3pgFUlUGxqfpHZlgYGyavtFJOgfG2X1ymGeOnqEgT7h8dS1Xrq23X6JmEVXlN7t6uGd7F++7YiXvfXkbdWVFfoc1p7kSvKc9PiJyNfBNIA+4RVW/4mV5xixW7/AETx89w9NH+tnaMcjhvlEK8gLUBAspK86nKlhAdbCQlpoS3rBxA8sqLKlnIxHh2gsauWxlNb/e0cVV//II562o4Mq1dVzYXMWGxgrqygozZr5+z2rwIpIHHADeAHQCzwLXqeqeePexGrxZDFVlYjrM2ekQM2F3JIRCWJ2v3WenQ4xPhRifnGF0coaxqRkGxqbpGZrgaP8ou08OMzI5w/rl5axbVs6ahjJaa4I2nNFwdirE7u4h9veMcLx/nOP9Y4gIq+pKWVVXSlttkOWVxdSWFlEVLKCypIDq0kKqg4XkBVLzIeBXDf6lwCFVPeIG8VPgbUDcBL8UHWfGGZuaQRVUIazKTFgJhcNMzSihsBJyP8wEZ4hUfkDIzwtQkCcERMgLCCIgOE9M9Id0Kp4qxYnd+V9RdX7argrT4TAzIWUmFCakStjdLyCQH3AeQ35egPyA81hEeO6xiDx/TiJ/p0NhpkPK1EyYyZkQUzNhpkJhQu6BAyIU5AUozA9QlO/8LcwPUBAIEAg8f1xwjhcOw0w4zEzYOfZMSJkJhwmH3e2qhNx9pt3HEXkM+QHn3OfnCaowE1YmZ0KcnQoxNhlifHqGSTeBn50KMTIx/VzzyemRKaaihrdFnifn8QsFAaGoII+i/AAlBXkUF+RRXBCgrDifiuICVlQGefk5dSyvKH7BAtczIWX47LSHz7bJFOsaylnXUA4479Ghs9N0DZ7l1PAEOzqHePTgaUYmphmZmGFscobhCWcKhKL8AM3VJayqK2NVXZAVVSXUlxdRWVJAsDCPwrw88gJCIABlRfk0VweTHruXCb4J6Ii63gm8LHYnEbkBuAGgtbV1UQX1Dk9w5VcfXtR9TeYrcRP4XKZmwkzNhBmypG08FBChqqSA6XCYsckQh/vGONw3ltB9t3/xjVSWJHdNAC8T/GyV3he1B6nqzcDN4DTRLKaghopijn3lzYu5qzHGZC0vx3R1Ai1R15uBLg/LM8YYE8XLBP8ssFZEVolIIfBO4B4PyzPGGBPFsyYaVZ0RkY8C9+EMk/y+qu72qjxjjDEv5Ok4MFW9F7jXyzKMMcbMzn5XbYwxWcoSvDHGZClL8MYYk6UswRtjTJZKq9kkRaQPOO53HAtUB5z2O4g0ZOclPjs3s7PzEt9c52alqtbPtiGtEnwmEpH2eBP95DI7L/HZuZmdnZf4FnturInGGGOylCV4Y4zJUpbgl+5mvwNIU3Ze4rNzMzs7L/Et6txYG7wxxmQpq8EbY0yWsgRvjDFZyhJ8AkTkahHZLyKHROSzc+z3P0VERSRnhnolcm5E5M9EZI+I7BaR21Mdox/mOy8i0ioiD4vIVhHZISLX+hFnqonI90WkV0R2xdkuIvKv7nnbISKXpjpGvyRwbt7lnpMdIvKEiFw070FV1S5zXHCmOj4MrAYKge3Axln2KwceBZ4CNvkdd7qcG2AtsBWodq83+B13mpyXm4G/dv/fCBzzO+4UnZtXAZcCu+Jsvxb4Dc6KcJcDT/sdcxqdm5dHvY+uSeTcWA1+fs8tHq6qU0Bk8fBY/wh8FZhIZXA+S+TcfBD4D1UdAFDV3hTH6IdEzosCFe7/leTIameq+ihwZo5d3gb8UB1PAVUi0pia6Pw137lR1Sci7yOcimTzfMe0BD+/2RYPb4reQUQuAVpU9depDCwNzHtugHXAOhF5XESeEpGrUxadfxI5L18C3i0inThrJnwsNaGlvUTOnYEP4HzTmZOnC35kiTkXDxeRAPAN4PpUBZRGEllYPR+nmeYqnBrHYyJyvqoOehybnxI5L9cBt6rq10TkCuBH7nkJex9eWkvk3OU0EXkNToJ/5Xz7Wg1+fvMtHl4OnA88IiLHcNoN78mRjtZEFlbvBH6pqtOqehTYj5Pws1ki5+UDwH8BqOqTQDHOhFK5LpFzl7NE5ELgFuBtqto/3/6W4Oc35+LhqjqkqnWq2qaqbThtY29V1XZ/wk2pRBZWvxt4DYCI1OE02RxJaZSpl8h5OQG8DkBENuAk+L6URpme7gHe646muRwYUtVuv4NKByLSCtwFvEdVDyRyH2uimYfGWTxcRP4BaFfV2Dduzkjw3NwHvFFE9gAh4FOJ1DwyWYLn5W+B74rIJ3GaIK5Xd3hENhORO3Ca6+rc/ocvAgUAqvptnP6Ia4FDwDjwfn8iTb0Ezs3fAbXAf4oIwIzOM8OkTVVgjDFZyppojDEmS1mCN8aYLGUJ3hhjspQleGOMyVKW4I0xJktZgjfGmCxlCd4YY7LU/weSyMSygAGDEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "NV =  3239 + 805\n",
    "NU =  3255 + 812\n",
    "RV = 1 / NV\n",
    "RU = (14 + 5) / NU\n",
    "VE = (RU - RV)/RU;\n",
    "\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV  + np.random.normal(0, 1, B) * np.sqrt(var_rv) + 1e-10\n",
    "RUs = RU  + np.random.normal(0, 1, B) * np.sqrt(var_ru) + 1e-10\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"Two-sided 95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.05))\n",
    "\n",
    "print(f\"One-sided 95% confidence interval of VE is [{CI_VE:.4f}, 1]\")\n",
    "\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try the Binomial parametric bootsrtap using the fact that outcomes are Bernoulli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "W-nfdQ9UzhzA",
    "outputId": "241c45ab-fbe6-46dc-ad5f-367bc69aa7e4",
    "papermill": {
     "duration": 0.539729,
     "end_time": "2021-01-04T05:11:01.934674",
     "exception": false,
     "start_time": "2021-01-04T05:11:01.394945",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall vaccine efficacy is 0.9471\n",
      "Two-sided 95% confidence interval of VE is [0.8000, 1.0000]\n",
      "One-sided 95% confidence interval of VE is [0.8333, 1]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEICAYAAABYoZ8gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxkVXnw8d9Te/XePd2zb+yLgKATxWgiCaJIFIy+iRKJYIzkjcaFNxpB/bgSNTFR9DVvFJWAGxoRiSYqsgiorAMzwAyzMcMsPTM93TPT3dVrrc/7x70lNTXV3VXVt+p2dz3fz6c+3XXvrXueW8tTp8499xxRVYwxxjSOgN8BGGOMqS9L/MYY02As8RtjTIOxxG+MMQ3GEr8xxjQYS/zGGNNgLPE3OBEZFZET/Y5jJiKyVkRURELu/Z+LyJUe7fsPRGRbwf3dIvIqL/bt7m+ziFzg1f7KLFNE5D9EZFBEHq1n2bNR/FqY2rDEX0Micp/7wYv6HctUVLVFVXd5vV8RuUpEfuP1fvNU9bWqeksZcaiInDzDvn6tqqd5EZeI3Cwi1xft/wWqep8X+6/AK4CLgJWq+pLCFSLyMhEZE5HW4geJyAYR+buCL9rRotubvQyy+PXx8rUoUdZxr02jssRfIyKyFvgDQIFLa1hOqFb7bgQL+PlbA+xW1bHiFar6ENALvKlwuYicBZwJ3FqwuMOtHORvP6hl0KZOVNVuNbgBHwN+C3wB+O+idTcDXwXuAkaA+4E1BesVeC+wCzgMfB4IuOuucvf7ReAocD3OF/hHgT1AP/AtoN3d/s3uftrc+68F+oCegrJOLojr/wE/B0bdcpYCNwCDwFbgvII4rwV2usfwDPCn7vIzgEkg6+5nyF0eBf4F2Asccp+D+BTPX9Dd9rAb/7vdWEPu+vuAv3b/P9l9Dofd7X/gLn/AfcyYG8ebgQtwkt6H3Ofh2/llBWXvBq5zj2kQ+A8gVvD8/6YoVnVjuBpIAym3vJ8W7O9VBc/BDcAB93YDEHXX5WP7e/d1PAi8fZr32HLgJ+774Fngne7ydxQ9/58s8dgPA/cWLftn4Hb3/7WFz3cZ7/d24JtuzPtx3pfBal+fotfig8BT7nbfBJbgvEdHgLuBzoLtf+i+rsPu/l/gLp/qtVkO/AgYAJ4D3luwr5cA64EEzvv1C37nFc/yk98BLNSb+0F8F/Bi9w23pGDdze6b9g/dRPClwmTifhh+BXQBq4HtPJ/krgIywHuAEBAH/sot70SgBbgd+HbB/r7rlrkIJ9m8rqiswsR/2I05BtzrfhjehpOIrwd+VfDYP3M/OAH3QzsGLCuIszhB3oCTqLqAVuCnwGeneP7+N84XzSp3+18xdeK/FfiIG0cMeEWp43PvX+A+f//kPvdxSiebTQVl/xa4fprjKn4Ory9av5vnE/+ngIeBxUAP8CDw6aLYPgWEgUuAcQoSW9F+78f5oo4B5+IkrwunirPosatw3per3fsBnC+dN7j311JZ4r8D+BrQ7B7bo8DfVPn6FL8WD+Mk+xU4X4hPAOe5r9+9wMcLtv8rnPdW/gt2Y9Hn7vqC+wHgcZxKWgTn87MLeI27/iHgL93/W4Dz/c4rXt18D2Ah3nDaV9NAt3t/K3BNwfqbge8X3G/BqZ2tcu8rcHHB+ncB97j/XwXsLSrvHuBdBfdPc8vPJ8kOnFr208DXih5bnLS+XrDuPcCWgvtn49bepzjujcBlBXEWfpkJzhfDSQXLXgY8N8W+7gX+d8H9VzN14v8WcCNOe3bxfkollhRuDb5gWXGyKSz7EmBnqeOa4jmcLvHvBC4pWPcanCaZfBwTFCRbnER3XMLBSdxZoLVg2WeBm6eKs8Q+7gY+7P5/Ec6Xfti9v9Y9rqGi2xkl9rMESFLw6w24HLeSUMXrU/xavLXg/o+Afy96j94xxfF1uPvP//o95rUBXsrxn6XrgP9w/38A+CTu53gh3ayNvzauBH6pqofd+99zlxXal/9HVUdxfq4vL7UepwlnqnW46/YUbR/C+UCiqkM4P4HPAv51htgPFfw/UeJ+S/6OiLxNRDaKyJCIDLn7755ivz1AE/B4wfa/cJeXspzjn4Op/APOF8ujbg+av5pmW4ABVZ2cYZvpnv/ZKPVaFe77iKpmCu6PU/CcF+3nqKqOFO1rRQWx3ILzaw7gL4HvqWq6aJtuVe0ouG0psZ81OL9QDha8tl/DqflD5a9PsbLekyISFJHPichOEUngfGnA1O/JNcDyfMxu3B/G/dzgNJmdCmwVkcdE5HUVxj1nLdQTW74RkTjw50BQRPrcxVGgQ0ReqKpPustWFTymBadJ4UDBrlYBm93/Vxet06JiD+C8ifNW4zQZHHL3fy7OT+BbgS8DF1d1cAVEZA3wdeBC4CFVzYrIRpwPeKkYD+N8SF+gqvvLKOIgBc8RzjGVpKp9wDvduF4B3C0iD6jqs1M9pIzyi8vOP/9jOF9guOUtrXDf+ddqqte2XAeALhFpLUj+q3Ha18t1O/D/ROSPgDfi1LarsQ+nxt9d9KUFVPX6VOsvgMuAV+Ek/XacczRTvSf34fziPKXUzlR1B3C5iARwnp/bRGSRljhhPt9Yjd97b8D5CX4mTrvruTgnO3/N87UrgEtE5BUiEgE+DTyiqoW1zA+KSKeIrALeB0zXm+JW4BoROcH9EvkMzgm0jIjEgO/g1GTeDqwQkXd5cJzNOB+kAQAReTtOjT/vELDSPT5UNYfzRfFFEVnsPmaFiLxmiv3/J/BeEVkpIp04J5JLEpE/E5GV7t1BN65sQRzVXKfwbrfsLpznLv/8Pwm8QETOdZ/bTxQ9bqbybgU+KiI9ItKN0778nUqDc98rDwKfFZGYiJyDU0P9bgX7GANuwzl5vUdV11cah7ufg8AvgX8VkTYRCYjISSLySqjZ61NKK84X0BGcL+fPFK0vLutRICEiHxKRuPuL4SwR+T037itEpMd97w65j8myAFji996VOG2Ee1W1L38DvgK8taD74PeAj+M08bwYeGvRfv4L58TTRuB/cHozTOUmnN4pD+CcjJ3EafsEp923V1X/XVWTwBXA9SJSspZTLlV9BqfZ6CGcD9TZOCdB8+7FqdX2iUi+yetDOCehH3Z/it+Ncz6ilK8Dd+Ik2idwaqdT+T3gEREZxTl5/D5Vfc5d9wngFven/J9XcIjfw0lmu9zb9QCquh3n5OvdwA6g+FqFbwJnuuXdUWK/1+P0FHkK55zLE/l9V+FynLb4A8CPcU5y3lXhPm7B+QXyrSnWDxX14/8/U2z3NpwTpPmeULcBy9x1tXh9SvkWTnPXfjeOh4vWH/PaqGoWeD1O5ew5nF+l38D5pQDOL+PNbtxfAt5SRhPhvCDuSQxTRyJyM04y/ugU6xU4pQY/hY0xxmr8xhjTaCzxG2NMg7GmHmOMaTBW4zfGmAYzL/rxd3d369q1a/0Owxhj5pXHH3/8sKoed5HkvEj8a9euZf36qroYG2NMwxKRkle8W1OPMcY0GEv8xhjTYCzxG2NMg7HEb4wxDcYSvzHGNBhL/MYY02As8RtjTIOxxG+MMQ2mZolfRG4SkX4R2VS0/D0iss2dgu2fa1W+Mcb4adfAKF/45Ta/wyipljX+myma4s+d4u0y4BxVfQHwLzUs3xhjfPPY7qN8+d5n+cWmvpk3rrOaJX5VfQBndqlCfwt8zp0JClXtr1X5xhjjp97BCc5e0c5H73iaofGU3+Eco95t/KcCfyAij4jI/fm5LUsRkatFZL2IrB8YGKhjiMYYM3v7Bsd52YmLWN4RZ+O+oZkfUEf1TvwhoBM4H/gg8J8iIqU2VNUbVXWdqq7r6TlucDljjJnT9g9O0N0apbMpwuHRxq7x9wK3q+NRIAd01zkGY4ypuYPDk3Q3R2iNhTg8mvQ7nGPUO/HfAfwxgIicCkRwZrY3xpgFI5dT+hNJFrVEaYuFGRhpkMQvIrcCDwGniUiviLwDuAk40e3i+X3gSrW5H40xC8zh0STN0SCRUID2eJj+xKTfIR2jZhOxqOrlU6y6olZlGmPMXLB/aILuligA7fEwG/Y1SI3fGGMa1YGhyecTf1O44U/uGmPMgrd/aJzO5gjg1PiPNPjJXWOMWfB6BydY5Cb+tliYxGSGbG7unM60xG+MMR7rHXy+jT8YEFqjIY6OzZ3mHkv8xhjjsYNDEyxqifzufkdTeE715bfEb4wxHhuaSNMafb7TZHt8bvXlt8RvjDEeG5nM0FSU+K3Gb4wxC1Qup4ynMjSFg79b1hqzxG+MMQvWaCpDLBwkEHh+/Mm2WMiaeowxZqEamczQHDl2UIT2pjD9CUv8xhizICUm0jRHg8csa4+HGbCmHmOMWZhGJjM0FdX4W2NhhsbTPkV0PEv8xhjjocREmqbIsTX+eDjIaDLjU0THs8RvjDEeGkmWSPyRIGOW+I0xZmFKTGSIl6jxj6UaIPGLyE0i0u9OulK87gMioiJi0y4aYxaUkck0sfCxiT8aCpDOKJlszqeojlXLGv/NwMXFC0VkFXARsLeGZRtjjC+GJtLHndwVEbe5J+tTVMeqWeJX1QeAoyVWfRH4B2DujFFqjDEeSUxkjmvjB2iKBBlJzo2ePXVt4xeRS4H9qvpkGdteLSLrRWT9wMBAHaIzxpjZS0ykaZ4i8c+Vnj11S/wi0gR8BPhYOdur6o2quk5V1/X09NQ2OGOM8UhiMk08cvx05vFIiNHJBkv8wEnACcCTIrIbWAk8ISJL6xiDMcbU1JQ1/nCQkTlS4z/+a6lGVPVpYHH+vpv816nq4XrFYIwxtZYoceUuQCwSWPg1fhG5FXgIOE1EekXkHbUqyxhj5orRZIam6PE1/rl09W7NavyqevkM69fWqmxjjPHL6GTpXj2xcHDh1/iNMabRTKaz5FSJBI9PrdFQkJHJBuzOaYwxC9nIZIbmaAgROW5dUyRIwmr8xhizsIxMpmkp0b4PThu/1fiNMWaBmapHDzgjdM6V7pyW+I0xxiMjk8cPyZwXt5O7xhiz8Djj9Exd47fEb4wxC8xYMkMsXDqtzqV+/Jb4jTHGI2OpDNFw6aaehhykzRhjFrrxVJZoaIoaf2TuzMJlid8YYzwyOpkhFipd44+FgkyksuRy/k9FYonfGGM8MjpNG38gIERDQcbT/s/CZYnfGGM8Mpqcuo0foCk6N3r2WOI3xhiPjKWmbuqB/Ale/6/etcRvjDEema47J0BTOMSI1fiNMWbhGE9miU/T1BOPBBd24heRm0SkX0Q2FSz7vIhsFZGnROTHItJRq/KNMabexlPZadv4o6EA43OgS2cta/w3AxcXLbsLOEtVzwG2A9fVsHxjjKmr8dT0TT2xcJCx5ALu1aOqDwBHi5b9UlXzX3cP40y4bowxC8JYKkuswWv8M/kr4Oc+lm+MMZ6aSGWn7dXjJP4FXOOfjoh8BMgA351mm6tFZL2IrB8YGKhfcMYYU4VcTklmskSnaeqJhAKMzYHxeuqe+EXkSuB1wFtVdcprl1X1RlVdp6rrenp66hegMcZUYSLtjNMTKDHtYl4sHGRsDtT4Sw8cXSMicjHwIeCVqjpez7KNMaaWxlKZadv3wZlwfUFfuSsitwIPAaeJSK+IvAP4CtAK3CUiG0Xkq7Uq3xhj6mmmPvwAsfDcOLlbsxq/ql5eYvE3a1WeMcb4qZwafyw0N5p67MpdY4zxwFhy+q6cANFwg57cNcaYhWgslSmjqSfYuN05jTFmoRlPTt+VE/KJ32r8xhizIDjz7U6fUhv6Ai5jjFloxpPTj8UPTo1/whK/McYsDGOpLJEpJlrPi4UDNvWiMcYsFGPJDNEZavyRYIBMNkfW5wnXLfEbY4wHRpMZ4jO08Yu4E677fILXEr8xxnhgponW8+IR/7t0WuI3xhgPjCenH5I5Lz4HLuKyxG+MMR4YnWGi9by5cBGXJX5jjPHAeBlj9QBuG78lfmOMmffGZ5h2MS8WDjBmJ3eNMWb+cxL/zCk1Gg4y7vOE65b4jTHGA+OpmfvxA8RCVuM3xpgFodwafyQU8H3YhlrOwHWTiPSLyKaCZV0icpeI7HD/dtaqfGOMqRdVZTJdXnfO6AKv8d8MXFy07FrgHlU9BbjHvW+MMfPaZDpHOBggEJh6ovW8aCi4cPvxq+oDwNGixZcBt7j/3wK8oVblG2NMvZQzCUteLLyAE/8UlqjqQQD37+KpNhSRq0VkvYisHxgYqFuAxhhTqYlUllikvMTvTL+4QNv4Z0tVb1TVdaq6rqenx+9wjDFmSmOpDLEZhmTOcyZcb6wa/yERWQbg/u2vc/nGGOO5ciZaz2vEGv9PgCvd/68E/qvO5RtjGkg6m+ODP3ySx3YXn270VrnDNcACr/GLyK3AQ8BpItIrIu8APgdcJCI7gIvc+8YYUxP/du+zbNg3xNXfWs9t6/fVrBynxl9mU88cmH4xVKsdq+rlU6y6sFZlGmNM3ta+BDc/uJt//NOzOTqW5IZ7dvCmF69EZOYul5UaT2WIBMtL/NGQDctsjDE18Z2H9vCaFyylqznCST0tpLM5th8arUlZ5Q7QBm6N3+d5d8tK/CLyIxH5ExGxLwpjzJynqtyztZ8Xr3EGBxARXrymkzs399WkvPFUZsaJ1vNi4cC8GZb534G/AHaIyOdE5PQaxmSMMbOy/dAoqsrKzvjvlr14dSc/33SwJuWNJbNEy078/rfxlxWpqt6tqm8FXgTsBu4SkQdF5O0iEq5lgMYYU6l7thzi3NWdx7Tnn7a0jQNDkxwYmvC8vNFkeSNzAoQCggKpTM7zOMpVdtONiCwCrgL+GtgAfAnni+CumkRmjDFV+uUzhzh3Zccxy4IB4dQlLTzVO+x5eWNlTrsITrNTPBxk3McuneW28d8O/BpoAl6vqpeq6g9U9T1ASy0DNMaYSoynMmztS3DGsrbj1q3siLPloPeJv5KTuwDxcJAxH5t7yu3O+Q1V/VnhAhGJqmpSVdfVIC5jjKnKpv0JVnc1lzzZuqqrmU0HalTjL7OpByAWCTDuY5fOcpt6ri+x7CEvAzHGGC88uW+IE3uaS65bs6iJrQdHPC9zLJUhWmZTDzgneP3s2TNtjV9ElgIrgLiInAfkz5S04TT7GGPMnPL43kFO7C6d+Je2xTgylmRkMk1rzLt+KeMVjNUD/g/bMFNTz2twTuiuBL5QsHwE+HCNYjLGmKo9tW+IV5+xpOS6QEBY3dXEtr4R1q3t8qzMsVSm7O6c4Pbl93GgtmkTv6reAtwiIm9S1R/VKSZjjKnK4dEkickMS9pjU26zuquJLR4n/ol0hTX+8Byu8YvIFar6HWCtiPyf4vWq+oUSDzPGGF881TvEKYtbCEwzHs+KjiY27/f2BO9Ehb16oiF/r96d6bdJvqGsBWgtcTPGmDnjyX1DrFk0/enHlZ1xth/y9gTveKr8K3cBIj4P1DZTU8/X3L+frE84xhhTvc0HEpy5rH3abZa0xegd9O7q3WxOSWVyZY/VA86E63O5xg+AiPyziLSJSFhE7hGRwyJyRa2DM8aYSmw5OMLqGWr8i5ojDI2nmfRohMz8JCzTNS8Vi4YCjM6DfvyvVtUE8DqgFzgV+GC1hYrINSKyWUQ2icitIjL1mRhjjCnDaDLDkbEky9qmTyeBgNDTGqV3cNyTcidSWeJlTrSeFwsHfW3qKTfx5zu8XgLcqqpVz2MmIiuA9wLrVPUsIAi8pdr9GWMMwLa+EVZ1NhEIzFzzXtIWY88RbxL/WCpLvIKLt8Dpzulnjb/cIRt+KiJbgQngXSLSA0zOsty4iKRxLgQ7MIt9GWMMW/sSrOoq77rSntYIe496lPgrGJkzLxaaB4O0qeq1wMtwaulpYAy4rJoCVXU/8C/AXuAgMKyqvyzeTkSuFpH1IrJ+YGCgmqKMMQ1k8/7EMePvT6enJcbuw2OelFvpAG0A0XCQMR8v4Krk98kZwJtF5G3A/wJeXU2BItKJ86VxArAcaC51olhVb1TVdaq6rqenp5qijDEN5JmDCdaUWeNf3BZlt2dNPeUPyZw3p6/czRORbwMnARuBfLQKfKuKMl8FPKeqA+6+bwd+H/hOFfsyxhhUlR39I6zqOqms7Ze0xTxr6plIZYlWWOOPhYOMp+d+G/864ExVVQ/K3AucLyJNOOcMLgTWe7BfY0yD2j80QSwcLHvgtcWtUfYPTZDLaVkng6fjtPFXWOOfD/34gU3AUi8KVNVHgNuAJ4Cn3Rhu9GLfxpjGtP3QCKvLbOYBp8bdHAnSP5KcddlVJX6fJ1wvt8bfDTwjIo8Cv3umVPXSagpV1Y8DH6/mscYYU2xb3ygrOso7sZu31G3uWTrNgG7lGKvi5K7fE66Xm/g/UcsgjDFmNrYcTFSc+Be1RD2ZeH1kMl3R7FvgjNWTzGQ9aWqqRrndOe8HdgNh9//HcJpqjDHGd9v6Rsruw5/X2RTmwPDsE//oZKbiGn9AhEgowIRHw0ZUqtyxet6J0y7/NXfRCuCOWgVljDHlymRz7D4yVnGNv6s5wn4PBmsbmcxUPGQD5Cdc96dnT7lnJN4NvBxIAKjqDmBxrYIyxphy7Tk6TldzpOJa96KWqCejdI4mM8QrLBvcLp0+9eUvN/EnVTWVvyMiIZx+/MYY46vtVTTzgDNK50EvmnqSlV/ABf7OwlVutPeLyIdxxte5CPgh8NPahWWMMeXZ1jfC8ip65ixqidI3PJshxxyjyQxNVTT1xML+9eUvN/FfCwzg9Lv/G+BnwEdrFZQxxpRr88EEKzsrr/G3xUJMpnOz7lY5lqz85C44ffn9Gpq5rO6cqpoTkTuAO/JDLRhjzFywrW+Ei85YUvHjRISe1ggHhic4qael6vLHU9nq2vh9vHp32hq/OD4hIoeBrcA2ERkQkY/VJzxjjJnaRCpLX2KSZR3VXYTlRV/+0WR1vXr8HJN/pqae9+P05vk9VV2kql3AS4GXi8g1NY/OGGOmsaN/hBUdcUKByk+ugtOl8+BQ9e38qspkOlvxBVyQH5p5bib+twGXq+pz+QWqugu4wl1njDG+2XpwhFVljsFfSmdTZFY1/vFUlkgoUNXVt7FQgNHJuZn4w6p6uHih285f3jB4xhhTI1v6Eqyo4sRu3qLmCL2zSPxOj55yR745VjwcZGSOJv5UleuMMabmnjmQYHVX9TX+RS2zq/FXe/EWQDwSJDGZrrrs2Zgp8b9QRBIlbiPA2fUI0BhjprL9kDPBerUWNc+uL/9YlSd2AeKRkG81/ml/o6hqdUdkjDE11j8ySSandDVHqt5HZ3OEQyPVJ/7RyVnU+MPBOdurxxhj5qTNBxKc2N2MSPXDGjdHgmSyWnUCdpp6qkuj8UiDJX4R6RCR20Rkq4hsEZGX+RGHMWb+2rR/mDWLmme1DxGhexZDNzgTrc+ixj9HT+7WypeAX6jq6cALgS0+xWGMmaee6h2qaLrFqXQ1RziUqC7xVzMWf15DNfWISBvwh8A3AVQ1papD9Y7DGDO/bd6f4ITu2dX4wUn81db4R5PZiufbzYtH5u4FXLVwIs6Ab/8hIhtE5BsictyrJyJXi8h6EVk/MGDDAxljnjc8nmZwPD3r+XIBOuJh+qqt8SfTs6rxz/Vhmb0UAl4E/LuqngeM4Yz+eQxVvVFV16nqup6ennrHaIyZwzYfGObEnmYCszixm9cxi6t3q519CyAcFFQhman/QG1+JP5eoFdVH3Hv34bzRWCMMWXZdGDYk/Z9cK7enU3ir7bGLyI0RYOM+TALV90Tv6r2AftE5DR30YXAM/WOwxgzf23YO8TaWfboyetsjlTd1DOWzNBUZeIHaIqEfOnZU90gE7P3HuC7IhIBdgFv9ykOY8w89PieQV595lJP9uX06klW9djZ1PjBHa8nWf9hG3xJ/Kq6EVjnR9nGmPntwNAEqUyOJW1RT/bXEQ+TmEiTyuSIVNhDZzZDNgA0Rfzpy29X7hpj5pXH9wxy2tLWWV2xWygQEDqawvRXMXTDWKr6IRvAv549lviNMfPK+t1HZzVVYimLWqJVXcQ1lszOrqkn4s/QzJb4jTHzymO7j3LqklZP99nVFOFgFRdxzbbG79f0i5b4jTHzxngqw67DY55csVuooylc8dW7uZwynszSNIs2/ljI2viNMWZajzznNPNUehJ2Jp1V1PhHJjPEItVNu5gX82kWLkv8xph5476t/Zy9ot3z/S5qibC/wou4hifStMZmNwOtk/jr353TEr8xZt64b9sA56zs8Hy/Xc0R+oYqq/EPTaRoic6uR3w8EmTE2viNMaa0fUfHSUymWbPIm6EaCi2q4urd4Yk0zbNo3wdo8mlMfkv8xph54f7tA7xwZYcnA7MV62yKcHg0STanZT9meCJNsxc1fkv8xhhT2i829XH2Su/b9wFCwQCtsRBHRssfumFoPE1zdHY1fr8mY7HEb4yZ8/oTk2zcN8SL13TWrIzulmhFPXuGJ9I0hWdf47fEb4wxJdyxcT8vOaGTaGh2NezpdDVX1qVzaDxN02ybeqyN3xhjjqeq/HB9Ly8/qbum5XQ2hekbLr9L5+B4atZNPS2xEAnrzmmMMcfasG+IxGSa05e11bScjqYIBypq6knREpldjT8SDKAKk+n6TsZiid8YM6d94ZfbeN05y2vSm6dQV4UzcQ2Oz75Xj4jQGgsxNF7fWr9viV9Egu5k6//tVwzGmLntsd1H2dE/ygWn1n7e7UUVtvEnPOjOCU5zz9BEatb7qYSfNf73AVt8LN8YM4elszk+9dNnuOzcFYSCtU9VXc1RDlZQ409MZGiZZRs/QEu0QWr8IrIS+BPgG36Ub4yZ+z5/5zbCQeGVdajtgzNez0AFF3ElJj2q8TdK4gduAP4ByE21gYhcLSLrRWT9wMBA/SIzxvju+4/u5ccb9vM3rzyp5m37eeFggLZYuKwJWdLZHMl0blZj8ec1R0MML/SmHhF5HdCvqo9Pt52q3qiq61R1XU9Pfb7xjTH+SmVy/Mud27jh7h1cd/HptM1y9MtK9bRGyxqlMzGRpjkW9GT6x6ZIsO41fj8mW385cKmIXALEgDYR+Y6qXuFDLMaYOSCXU+7Z2s/nfr6F9niYj73+TDqbInWPo7slyv7BCX5v7fTbDYW02WAAAA/GSURBVE2kaY1686XUHA0xOF7fGn/dE7+qXgdcByAiFwAfsKRvTGMaTWb40eP7uPnBPQQELjt3BevWdHo2kXqluprDZdX4hyfSsx6SOa8lGmJwbOHX+I0xDa5veJKv3r+T25/o5QXL27jipas5Y1mbbwk/b1FLlL1HxmfczhmZ05vhI1qiIfYdnblML/ma+FX1PuA+P2MwxtTP0HiKG+7ezo+e2M8rT+3hH//0bLpbon6H9TvdLVEe7Ds843bDHly8ldfSCE09xpjGo6r8eMN+Pv3fz/CStV18/n+9kPZ4fU/clqOnpbyTu8MT6VlNsl7I6dVjTT3GmAVkeDzNB374JM8OjPCBV5/GiT0tfoc0pe6WKH2JSVR12mano2NJT2v89U78NlaPMaZmtvYluOTLvyYUFD556VlzOumDMz5+OBjg6Nj0TS8Hh5Oe9TpqiYZITNR3aGZL/MaYmrhvWz9vufFh3nDeCt72srWE6zDsghcWt0bpHZy+uacvMelZ4o+FA6SyOZKZ+o3QOT9eCWPMvHLHhl6u+cFGrnnVqbzi5NqOo++17jLa+fsTk3Q2eXOOQkRoi9W3uccSvzHGU7c+upfr/2cL1732DE5d0up3OBVb3Bpl95Gxabc5lJiks9m7C8xaYiGG63j1rp3cNcZ45ruP7OGGu3fw4UvOYFl73O9wqrKkPcaz/aNTrk9nc4xMZmj3cDiJlmiYwTomfqvxG2M88YPH9nLD3Tv4yDxO+gDL2uPsGpi6xj8wkqSjKUwg4N3FZi3RIEN17Mtvid8YM2u3Pb6Pz9+5jeteezpL2mJ+hzMrS9ti7JmmqedQYpIuD5t5wOnLX8+B2izxG2Nm5Y4NvXz2Z1u59rXzu6af19kUZiKdnfJk66GEd10581qjIY7M0IXUS5b4jTFV+/ETvXzqv7fwoYtPZ0XH/E/64PSyWdERZ/fh0rX+Q4lJOjy+6rijKULfcPmzf82WJX5jTFX+8zGn9861F5/Oqq4mv8Px1LL2GM9Nkfj7hidp8zjxdzaFK5rvd7asV48xpiKqytcf2MU3f/uccyJ3gdT0Cy1ui7FroHTPnoOJCZa0enseo7M5UtbMX16xxG+MKVsmm+PjP9nMb3Yc5qN/cuacGlnTS0vbYjw7ReI/NJzk9CVtnpbX2RShfyTp6T6nY4nfGFOWI6NJ3vXdJ0hlc3zs9WfSFFm46WNZe5x7t/aXXNc/4u3FW+Ak/sOjyRkHh/OKH3PurhKRX4nIFhHZLCLvq3cMxpjKPLzrCJd8+dcsa4/zgYtOW9BJH2BFR5w9R8ZJZ3PHresfSXo2XENeJBQgHg7OODicV/x49TLA36vqEyLSCjwuInep6jM+xGKMmUYyk+WLd23nh+t7+es/OJFzV3X4HVJdxCNBFrdF2dY3wlkr2n+3fHAsRS6nnk27WKirOcKhRJJFdWg+q3uNX1UPquoT7v8jwBZgRb3jMMZM7+neYV735d/wxN4hrn/DWQ2T9PNO6mlh476hY5Y9vX+YE3qaa9Ic09Uc4dBIfU7w+tqdU0TWAucBj/gZhzHmealMjs/fuY2//OYjXHTmEt5/4Sl0eHzB0nxwQnczT+wZPGbZpv3DrF3UXJPyOpoi9NepZ49vDXUi0gL8CHi/qiZKrL8auBpg9erVdY7OmMa0+cAw7//+RjqawnzmjWd7foXqfHJSTwvf2LrrmGVP9g5xyuLajDjaHg9zKFGfnj2+1PhFJIyT9L+rqreX2kZVb1TVdaq6rqenp74BGtNgsjnlK/fu4K1ff4QLz1jCNa86taGTPsCqrjgHhycZmXx+6Ian9ydY212bGn9nU5gDZcz364W61/jFaRz7JrBFVb9Q7/KNMcfad3Sc931/A+ms8uk3nLVg++ZXKhQIcEJ3M0/3DvP7J3czNJ5iaDzFsvbaDELX2RRh1+HBmTf0gB81/pcDfwn8sYhsdG+X+BCHMQ3vvzbu5/X/9zecsayNa197uiX9IqcvbeUXm/sA2LQ/wQndzQRq1M++nlfv1r3Gr6q/AWp/hYIxZkqJyTQf/fEmNuwd5B8uPp0TatR8Md9deMYSrr39Kd7/qlN59LkjNTuxC+7Vuwu5jd8Y45+Hdx3hNV98gMl0lk+/4SxL+tPoao7w0hO6uOYHG/n2w3u44LTanW/siIcZS2UYTWZqVkaeJX5jGsR4KsMnfrKJd3/3Ca44fw1vf/kJRENBv8Oa8y45axlP7hvimledysrO2o1CGgg4w0HvnGbaR68s7OuujTEAPLB9gOtuf5qTepr57BvPptXD+WIXumUdcf7trS+qWdt+oRUdcXb0j/LCGl8sZ4nfmAVs/9AEn/rpZp7cN8yVv7+Gc1d1+h3SvFSPpA+wvCPOtr6Rmpdjid+YBWg8leGr9+3klof2cNEZS/inN51DJGQtu3Pdio44j+2pfZdOS/zGLCC5nPKjJ3r5/J3bOHVJK5++7Cx6Wq2L5nyxojPOrY/tq3k5lviNWSAe2nmET/50MwB/90cnc8qS2gwtYGpncWuMI6NJxlOZmg59bYnfmHlu58Aon/mfLWw6MMyb163i/BMX1WUyD+O9YEBY3hFn18DYMcNBe80SvzHzVO/gOF++Zwe/2NzH685ezttettba8ReAFR1xth8ascRvjHGoKo/vGeSWh3Zz/7YB/vj0Jfzrn51bk4lBjD9WdzXx+J5B3viilTUrw94txsxhqsrASJKN+4Z4cOdhfrn5EIGAcMFpPXzhz8+l2RL+gvPCVR186e7tNZ1/1941xvhsIpXl2f5Rdh8Z48DQBAeGJjg4PEnv0AS9R8cBOLGnhdOWtPLeC09hdVeTteEvYKs646Rzys6BMU5e3FKTMizxG1Nng2MpfrvzML/ZcZj1uwfZOzjOio44S9tjLGqO0BGPcPrSVl5+cjeLW6O0x8OW6BuIiHDuqg5+tfWQJX5j5quJVJYN+wZ58Nkj3L99gJ0Do5y5rI0zlrVx1cvXsqariVDQTsqa552zsp27t/Tzzj88qSb7t8RvjEeyOWX/4AQ7D4+ys3+ULQcTbNqfYM/RMdZ0NXP6slYufeFyTlvaStgSvZnGWcvb+dr9u+gdHK/JwHCW+I2pUDanPHd4lC0HR9h6MMG2QyPsHBhj/+AEHU1hlrvNNsvb41xx/mpWdzVbN0tTkVg4yMUvWMLnfr6Vr/zFizzfvy+JX0QuBr4EBIFvqOrn/IjDmGK5nDKezjKWzJCYSHN0LEVfYpLewQl29I+w/dAouwZG6WqKsGZRM8s7Ypy5rJ2LzlzK0raYJXjjmT85Zzkf+OGTbNg7yHmrvR1cz485d4PAvwEXAb3AYyLyE1V9pt6xmPknl1PSuRzZnKL6/PKACCLP/wWnZp7O5khmcownswxNpDg6lqJ/JEl/YpL9bu+ZgZEkR8dSJCbTjCezRMMB4uEgLbEQrdEwnc1hupoiLG2Pc96qTlZ2xmt6Ob0x4NT6X3laD3du7pv/iR94CfCsqu4CEJHvA5cBNUn8wxNpDiUmUQXFyRSqzi2nTvLIqpLNKbmCv7mcs31xcgkEIChCICAEA+L87yabwsQj7uySxZ0xZts3Q0stc4+tMOZsTn93XOlsjnRWSWdypLM5Uu79bC5HJqfk1OkvXirOnLuvTC5HKuPcku4tlcmScRNwICAEBcLBAOFggFAwQL4ZO5uDVCbHRDrDWDLDSDLD6KQz09B4KstkOksq48bixp3LPf+6eOH5RB6ivSlMZ1OEzqYwJy9uYd2aLtrjIZojIeKR4IxD8GaySmIi7UlcxkwnlwP15iNwDD8S/wqgcPi5XuClxRuJyNXA1QCrV6+uurC//c7jPLjzSNWPN/4LihAKCuFggGBAqv7yTGdyHM04tf7nGPM0RmNq5ZqLTvF8n34k/lKf2+O+01T1RuBGgHXr1lX9nfe9d55f7UONMWZB8uNMVC+wquD+SuCAD3EYY0xD8iPxPwacIiIniEgEeAvwEx/iMMaYhlT3ph5VzYjI3wF34nTnvElVN9c7DmOMaVS+9ElT1Z8BP/OjbGOMaXR2tYkxxjQYS/zGGNNgLPEbY0yDscRvjDENRoov1Z+LRGQA2FOwqBs47FM4frDjXdga6Xgb6VjB/+Ndo6o9xQvnReIvJiLrVXWd33HUix3vwtZIx9tIxwpz93itqccYYxqMJX5jjGkw8zXx3+h3AHVmx7uwNdLxNtKxwhw93nnZxm+MMaZ687XGb4wxpkqW+I0xpsHM6cQvIheLyDYReVZEri2x/ioRGRCRje7tr/2I0yszHa+7zZ+LyDMisllEvlfvGL1Uxuv7xYLXdruIDPkRpxfKONbVIvIrEdkgIk+JyCV+xOmVMo53jYjc4x7rfSKy0o84vSAiN4lIv4hsmmK9iMiX3efiKRF5Ub1jPI6qzskbzpDNO4ETgQjwJHBm0TZXAV/xO9Y6Hu8pwAag072/2O+4a3m8Rdu/B2cIb99jr9FreyPwt+7/ZwK7/Y67xsf7Q+BK9/8/Br7td9yzON4/BF4EbJpi/SXAz3FmHzwfeMTvmOdyjf93k7KragrIT8q+UJVzvO8E/k1VBwFUtb/OMXqp0tf3cuDWukTmvXKOVYE29/925vesdOUc75nAPe7/vyqxft5Q1QeAo9NschnwLXU8DHSIyLL6RFfaXE78pSZlX1Fiuze5P59uE5FVJdbPF+Uc76nAqSLyWxF5WEQurlt03iv39UVE1gAnAPfWIa5aKOdYPwFcISK9OHNVvKc+odVEOcf7JPAm9/8/BVpFZFEdYvND2e/1epnLib+cSdl/CqxV1XOAu4Fbah5V7ZRzvCGc5p4LcGrA3xCRjhrHVSvlHG/eW4DbVDVbw3hqqZxjvRy4WVVX4jQNfFtE5vLnczrlHO8HgFeKyAbglcB+IFPrwHxSyXu9LubyG2vGSdlV9YiqJt27XwdeXKfYaqGcSeh7gf9S1bSqPgdsw/kimI/KOd68tzB/m3mgvGN9B/CfAKr6EBDDGeBrPirns3tAVd+oqucBH3GXDdcvxLqq5L1eF3M58c84KXtRO9mlwJY6xue1ciahvwP4IwAR6cZp+tlV1yi9U87xIiKnAZ3AQ3WOz0vlHOte4EIAETkDJ/EP1DVK75Tz2e0u+EVzHXBTnWOsp58Ab3N795wPDKvqQT8D8mXO3XLoFJOyi8ingPWq+hPgvSJyKc5PxKM4vXzmpTKP907g1SLyDJAFPqiqR/yLunplHi84TSDfV7d7xHxU5rH+PfB1EbkGpxngqvl6zGUe7wXAZ0VEgQeAd/sW8CyJyK04x9PtnqP5OBAGUNWv4pyzuQR4FhgH3u5PpM+zIRuMMabBzOWmHmOMMTVgid8YYxqMJX5jjGkwlviNMabBWOI3xpgGY4nfGGMajCV+Y4xpMP8fQmzCmFjcjZAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "NV =  3239 + 805\n",
    "NU =  3255 + 812\n",
    "RV = 1 / NV\n",
    "RU = (14 + 5) / NU\n",
    "VE = (RU - RV)/RU;\n",
    "\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = np.random.binomial(NV, RV, size=B) + 1e-10\n",
    "RUs = np.random.binomial(NU, RU, size=B) + 1e-10\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"Two-sided 95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.05))\n",
    "\n",
    "print(f\"One-sided 95% confidence interval of VE is [{CI_VE:.4f}, 1]\")\n",
    "\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More Accurate Confidence Intervals\n",
    "\n",
    "The table from the FDA uses exact inference by inverting tests based on the exact binomial nature of the outcome variable. This method is known as the Cornfield Procedure to find the exact confidence interval on the estimate of vaccine efficacy.\n",
    "\n",
    "Another typical more accurate approximation is to approximate the log of the risks by a normal distribution, construct confidence intervals for the quantity $\\log(RV) - \\log(RU)$ and then invert. This is for instance method C in [this paper](https://www.jstor.org/stable/pdf/2530610.pdf) and is implemented in the statsmodels package (and also in the `scipy.stats.contengency_table` package)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>        <th>Estimate</th>  <th>SE</th>     <th>LCB</th>    <th>UCB</th>  <th>p-value</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Odds ratio</th>        <td>0.053</td>      <td></td>  <td>0.027</td>  <td>0.104</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Log odds ratio</th>   <td>-2.930</td> <td>0.342</td> <td>-3.601</td> <td>-2.260</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Risk ratio</th>        <td>0.054</td>      <td></td>  <td>0.028</td>  <td>0.105</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Log risk ratio</th>   <td>-2.922</td> <td>0.342</td> <td>-3.593</td> <td>-2.252</td>   <td>0.000</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tb = sm.stats.Table2x2([[9, 19965 - 9], [169, 20172 - 169]])\n",
    "tb.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall vaccine efficacy is 0.9462\n",
      "Approximate 95% confidence interval of VE: [0.8948, 0.9725]\n"
     ]
    }
   ],
   "source": [
    "print(f\"Overall vaccine efficacy is {1 - tb.riskratio:.4f}\")\n",
    "lb, ub = tb.riskratio_confint(alpha=.05)\n",
    "print(f\"Approximate 95% confidence interval of VE: [{1 - ub:.4f}, {1 - lb:.4f}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that this matches much more accurately the table from the FDA."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "[tbd]vaccine-rcts.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  },
  "papermill": {
   "duration": 8.643307,
   "end_time": "2021-01-04T05:11:02.054671",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-01-04T05:10:53.411364",
   "version": "2.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
